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Abstract 



Low Reynolds-number, mildly curved, turbulent channel flow has been sim- 
ulated numerically without subgrid scale models. A new spectral numerical method 
developed for this problem was used, and the computations were performed with 
2 million degrees of freedom. A variety of statistical .and structural information 
has been extracted from the computed flow fields. These include mean velocity, 
turbulence stresses, velocity skewness and flatness factors, space-time correlations 
and spectra, all the terms in the Rcynolds-stress balance equations, and contour 
and vector plots of instantaneous velocity fields. 

The effects of curvature on this flow were determined by comparing the 
concave and convex sides of the channel. The observed effects are consistent with 
experimental observations for mild curvature. The most significant d'.' rencc in 
the turbulence statistics between the concave and convex sides was in the Reynolds 
shear stress. This was accompanied by significant differences in the terms of the 
Reynolds shear stress balance equations. In addition, it was found that stationary 
Taylor-Gortler vortices were present and that they had a significant effect on the flow 
by contributing to the mean Reynolds shear stress, and by affecting the underlying 
turbulence. 

Turbulence statistics were found to be in qualitative agreement with the 
large-eddy simulation of a plane channel performed by Moin & Kim (J. Fluid Mech. 
118, 341, 1982). Near the walls, the flow consists of alternating high and low speed 
streaks, with mean spanwise spacing of 100 wall units. It was also found, however, 
that near the wall, the velocity fluctuations normal to the wall are dominated by 
small, intense regions that are not significantly elongated in the streamwise direc- 
tion. Streamwise vortices were observed near the walls and were found to occur 
most often as a single vortex, rather than in vortex pairs. 
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1. Introduction 


Turbulent flow over curved walls is of considerable engineering interest. It occurs, 
for example, in turbomachinery and on airplane wings and in many other applica- 
tions. However, current methods for predicting these flows are quite inadequate, as 
is evidenced by their poor performance when applied to the relatively simple cur- 
vature cases in the 1980-81 AFOSR-HTTM- Stanford Conference on Complex Tur- 
bulent Flows (Kline, Cantwell, & Lilly 1982). One of the reason for this difficulty' 
is what Bradshaw (1973) calls “the surprisingly large effect exerted on shear-flow 
turbulence by curvature of the streamlines in the plane of the mean shear.” He 
notes that curvature effects are often an order of magnitude greater than would be 
predicted by using dimensional arguments. This poor understanding of the effects 
of curvature greatly hinders modeling efforts. 

1.1 Curvature Effects 

The effects of curvature on fluid flow have been under experimental and theo- 
retical investigation for some time. An inviscid stability analysis first performed 
by Rayleigh (1917) indicates that flow over a concave curved surface is unstable 
and, t /ersely, that convex curved flows are stable. Go* tier (1940) performed a 
viscous stability analysis showing that laminar flow over a concave surface is un- 
stable at suflicientlv high Reynolds number. This instability leads to a system of 
large longitudinal roll cells. These so-called Taylor- G or tier cells were later observed 
experimentally (Gregory & Walker 1950). 

Early experimental studies of the effect of curvature on turbulence (Wilcken 1930, 
Wattendorf 1935) revealed changes in mean-flow properties much larger than had 
been predicted by mixing-length arguments. Boundary layers wer* observed to grow 
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much faster on concave surfaces than on flat ones and, conversely, to grow slower on 
convex surfaces (Wilcken). Wall shear stresses were also greatly affected, increasing 
on a concave wall and decreasing on a convex wall. In a fully developed curved 
channel with strong curvature {6/R w 0.1) Wattendorf observed a constant angular 
momentum ( rU ) profile over a large portion of the central flow region; however, 
Hunt Sc Joubert (1979) found that for a weakly cui-ved, fully developed channel 
( 6/R « 0.01) this constant angular momentum region was not present. 

Sufficiently close to the curved walls, mean velocity profiles have been observed 
to obey the “law of the wall.” This has been the case for both concave and convex 
curved flows in boundary L. ers (So Sc Mcllor 1973, 1975, and others) and in fully 
deveh ped curved channel flow (Ellis Sc Joubert 1974). At larger distances from the 
wall, the mean velocity of a convex wall layer exceeds that of the flat-wall profile and 
the mean-velocity profile of a concave wall layer lies bellow the flat wall profile when 
plotted in law of the wall coordinates. The point at which these deviations occur and 
their magnitudes are dependent on the curvature parameter {6/R). For sufficiently 
weak curvature {6/R ^ 0.01), these deviations occur beyond the logarithmic region 
(Hunt Sc Joubert 1979; Hoffman Sc Bradshaw 1978); as an example, mean velocity 
profiles from Hunt and Joubert are shown in Figure 1.1. It has been suggested 
(Hoffman Sc Bradshaw 1978) that the flat-plate law of the wall applies where y/R 
is small (y is distance from the wall). 

Turbulence quantities are also affected by curvature, as would be expected from 
Rayleigh’s analysis. In the strongly curved convex boundary layers {6/R « 0.05 - 
0.1) of So Sc Mellor (1973) and Gillis Sc Johnston (1983) the turbulence virtually 
vanished in the outer half of the boundary layer. The point at which the turbulent 
shear stress fell to zero was well within the velocity gradient layer. This led Gillis 
and Johnston to propose that outer layer scaling for these flows should be based 
on the thickness of the shear-stress layer instead of on the classical boundary-layer 
thickness. Gillis and Johnston also observed that when shear stress, normalized by 
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friction velocity (- uv/u 2 ), was plotted against y/R all of their data, as well as the 
data of So and Mcllor, collapsed to a single curve. 

Boundary layers on strongly curved concave surfaces ( 6/R a 0.1), as studied by 
So and Mcllor (1975), show dramatic increases in turbulence activity. In the outer 
layer, intensities and shear stress were observed to be about twice as large as in 
flat wall boundary layers at similar conditions. The point at which the shear stress 
fell to zero was at y/6 a 1.1. However, both turbulent shear stress normalized 
by turbulence energy (— uv/q 2 ) and energy normalized by friction velocity (q 2 /u 2 ) 
showed reasonable agreement with flat-wall data near the wall. This is further 
evidence of the similarity of the near-wall regions. 

In the weakly curved boundary-layer cases {6/R a 0.01) of Hoffman & Bradshaw 
(1978) and Muck (1982), smaller changes in turbulence quantities were observed. 
Changes in turbulence intensities were 10% to 20%, increasing on the concave wall 
and decreasing on the convex wall. Turbulent shear stress increased or decreased 
about 10% relative to the flat-wall case, with most of the change occurring in the 
outer layer. Even these modest changes are noteworthy since an order of magnitude 
analysis of the Reynolds-stress transport equations predicts changes an order of 
magnitude smaller for this mild curvature. In these investigations changes in third 
and fourth order statistics of order one were also observed. 

In fully developed curved channels, similar changes in turbulence quantities have 
been observed (Eskinazi & Yeh 1956; Hunt h Joubert 1979). Turbulent intensities 
are larger on the concave side and smaller on the convex side. Also, the point where 
the turbulent shear stress is zero is displaced significantly toward the convex wall, 
and the wall shear stress is larger on the concave side than on the convex side. 

Tani (1962) has suggested that there is a turbulent analog to the laminar Taylor- 
Gortler vortices. He was led to this proposal after observing stationary spanwise 
variations in mean velocity in a concavely curved boundary layer. Similar obser- 
vations have since been made by many researchers (Patel 1968; So & Mellor 1975; 
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Meroney Sc Bradshaw 1975; and others). Evidence of turbulent Taylor-Gortler ceils 

has also been found in fully developed channels (Hunt Sc Joubert 1979). These 

longitudinal vortices give rise to spanwise variations in boundary- layer thickness 

and skin friction. Boundary-layer thickness is greatest at the boundaries between 

the assumed vortices where the motion is away from the wall (outflow), and skin 

friction is lowest there. Turbulence quantities in the outer layer are also affected by 

these large longitudinal structures. In general, the effects of concave curvature on • 

turbulence quantities are greater at the outflow boundaries between the postulated 

roll cells (So Sc Mellor 1975; Muck 1982). 

Many researchers have observed a repeatable stationary pattern of spanwise vari- 
ations. This repeatability has been attributed to upstream disturbances (Meroney 
Si Bradshaw 1975). In an attempt to impose a more regular pattern of variations, 

Muck (1982) placed regularly spaced vortex generators upstream of the curved 
section. The resulting weak longitudinal vortices in the upstream boundary layer 
were amplified by the curvature, serving to “lock in” the positions of the turbulent 
Taylor-Gortler cells. 

Interestingly, Jeans Si Johnston (1982) did not observe a stationary pattern of 
roll cells in their flow visualisation study on concave curvature, presumably because 
of a lack of persistent upstream disturbances. They observed large-scale roll-like 
structures (they referred to them as sweeps and ejections) which appeared randcmly 
in time. These structures drifted in the spanwise direction and had streamwise 
extent as small as “several boundary- layer thicknesses.” Barlow (1983), using the 
same experimental facility, observed that the appearance of these structures was not 
entirely random, and that their extent and persistence were dependent on upstream 
conditions, for example, the condition of screens in the water channel. Barlow also 

Jl 

placed vortex generators upstream of the curved section and was able to make 
the roll-cell pattern stationary, with rolls extending the entire length of the curved 
section. 
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To model the effects of ..tuvaiure on the mean-flow properties, Bradshaw (1973) 
suggested a correction factor for the apparent mixing length in analogy to the 
Monin-Obouhkov formula for buoyant Bows: 

<„ ' au/a, ’ 

where l and Iq are the corrected and uncorrected mixing lengths, and (3 is an 
empirical constant of the order o f 10. This model has met with limited success in 
cases of weak curvature. More complicated schemes involving the solution of the 
modeled Reynolds stress transport equations have also been used (Irwin Sc Smith 
1975; Gibson,, Jones, Sc Younis 1981). The Reynolds-stress transport equations have 
additional production, convection, and diffusion terms arising from the curvilinear 
coordinate system. These terms do not appear in scalar equations, such as those 
used in models based on the k-e equations. In addition, it has been suggested 
that curvature terms which arise naturally in the model of Launder, Reece, Sc Rodi 
(1974) for the pressure strain correlation may account for observed curvature effects 
(Launder et al.). Thus, by solving the Reynolds-stress equations, the presence c r 
curvature is reflected in the equations being solved, rather than explicitly added to 
the turbulence models being used. This approach has enjoyed reasonable success. 

Modeling efurts for curved flows have been hindered because the mechanism by 
which curvature induces the dramatic changes noted above is not well understood. 
The study reported here was undertaken in an attempt to improve this understand- 
ing by using numerical simulation to provide data that are not normally available 
from experiments. Numerical simulation of a turbulent flow can provide the turbu- 
lent velocity field as a function of space and time, which can be used to compute any 
quantity of interest. For example, all of the terms in the Reynolds-stress transport 
equation can be computed, as can two-point space-time correlation functions. In 
addition, the instantaneous velocity and pressure fields can be examined to gain 
information about the turbulent structures they contain. 
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1.2 Turbulent Flow Simulation 


The numerical simulation of turbulent flows is a technique for investigating turbu- 
ence, which has recently become feasible because of the development of large scale 
computers. These simulations are proving to be valuable supplements to labora- 
tory measurements. In one type of simulation (direct simulation), the Navier-Stokes 
equations (usually incompressible) are numerically solved for a turbulent flow field. 
However, such a simulation must resolve all important spatial and temporal scales 
in the flow. One of the characteristics of turbulence is that the range of scales 
in the flow increases rapidly with Reynolds number. Thus, a direct simulation is 
necessarily limited to low Reynolds numbers. 

Another method, the large-eddy simulation (LES), alleviates the low-Reynolds- 
niunber restriction at the expense of including some modeling. In LES, the large- 
scale motions are explicitly computed, but the effect of the unconputed small scales 
is modeled. This strategy relies on the fact that the large-scale structures vary 
significantly from flow to flow, whereas the small-scale eddies are thought to be 
universal (Chapman 1979). This universality of the small scales makes them much 
more amer oie to modeling than the turbulent flow as a whole. However, there 
ire difficulties with the LES approach when computing the near-wall region of 
a wall-bounded flow. In this region, the important large eddies (the streaks of 
Kline, Schraub, & Runstadler 1967) are in fact quite small relative to the flow as a 
whole, and they become smaller as the Reynolds number increases. Thus, the LES 
technique for wall-bounded flows is also limited to low Reynolds numbers when the 
viscou sublayer is computed. However, this limitation can be alleviated somewhat 
b” the use of fine grids embedded near the walls. For more details on the techniques 
of turbulent flow simulation, see the review article by Rogallo & Moin (1984). 

3oth direct simulations and LES put a great burden on computing hardware, and 
because ol the great computational efforts involved they are limited to very simple 



flow situations. To date, flows that are at most one dimensional in the mean have 
been computed. The other two dimensions are assumed to be homogeneous, and 
periodic boundary conditions are used. The use of periodic boundary conditions 
greatly simplifies the computation. The current limitations of these simulations 
to one dimensional, low-Rcynolds- number turbulent flows makes them impractical 
for engineering calculations. However, they are being used as a research tool in 
fundamental studies of turbulence. 

Unbounded turbulent flows have been computed quite successfully using both 
direct simulation and LES. Homogeneous isotropic turbulence has been computed 
by Orszag Sc Patterson (1972) and by many others using direct simulation, and by 
Kwak, Reynolds, Sc Ferziger (1975), Shannan, Ferziger, Sc Reynolds (1975), Man- 
sour et al. (1979), and Antonopoulos-Domis (1981) using large eddy simulation. 
Rogallo (1981), Fciereisen, Reynolds Sc Ferziger (1981), and Shirani, Ferziger, Sc 
Reynolds (1981) have performed direct simulations of homogeneous turbulence un- 
der the influence of various mean strains. Large eddy simulations of these flows 
have also been performed (Kwak, Reynolds, Sc Ferziger 1975; Shaanan, Ferziger, 
Sc Reynolds 1975; and Bardina, Ferziger, Sc Reynolds 1983). The results of these 
computations have been Uoed as test cases for turbulence models. In addition, 
simulations have been performed for time-developing unbounded plane shear layers 
(Mansour et al. 1978; Cain et al. 1981; and Riley Sc Metcalf 1980). 

Simulation of wall-bounded flows has proved to be more difficult, because the 
energy-producing turbulence structures near the wall are small compared to the 
overall dimensions of the flow. Deardroff (1970) and Schumann (1973) have per- 
formed simulations of fully developed plane channel flow using LES. They were 
able to compute the the flow far from the wall, but the effects of the walls and the 
near- wall turbulence were modeled. Moin & Kim (1982), using LES, actually com- 
puted the near-wall flow, but their numerical resolution was insufficient to properly 
resolve the near-wall turbulent eddies. However, their calculations did reproduce 
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the experimentally observed structure of the flow near the wall, though at scales 
larger than reported experimentally. A near- wall subgrid scale model was included 
to account for the missing production of the under-resolved eddies. The simula- 
tions being reported here were performed using direct simulation, at sufficiently 
low Reynolds number for all important turbulence scales to be resolved. 


1.3 Numerical Methods 

Turbulent flow simulations are often performed using spectral methods. Spectral 
methods are used, because for sufficiently smooth fields, they have a very high 
formal accuracy. This is particularly important in three-dimensional simulations, 
in which the number of modes that can be used in each spatial direction is severely 
limited. In addition, the periodic boundary conditions often used in these problems 
make spectral methods based on Fourier expansions quite natural and easy to apply. 
However, the application of spectral methods to wall-bounded flows is considerably 
more difficult. Fourier expansions are not appropriate for the direction normal to 
the wall because the imposition of the no-slip boundary condition severly degrades 
their accuracy. There is also a fundamental numerical difficulty associated with 
the incompressible Navier-Stokes equations in wall-bounded domains when solved 
using spectral methods. The difficulty stems from the continuity equation and 
the no-slip boundary conditions, which appear as constraints to the Navier-Stokes 
equations. When the dynamic equations are time-differenced, the continuity and 
boundary constraints must be imposed on the velocity field at each time-step. Moin 
& Kim (1980) have shown that when spectral methods are used, the most common 
explicit time-advancement scheme leads to meaningless calculations, because the 
continuity and boundary conditions cannot be properly enforced. They suggest 
implicit time- differencing of the viscous and pressure terms to allow the imposition 
of the constraints. 
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Several schemes have been developed for solving the Navier- Stokes equations in 
wall-bounded domains, asing spectral methods. Most of these computations have 
used Fourier expansions in two space dimensions. In the method described by 
Moin & Kim (1980), the velocity and pressure arc expressed in terms of Chcbyshev 
polynomials (and Fourier functions). The momentum equations are time-differenced 
with the viscous and pressure terms treated implicitly. The resulting equations are 
solved simultaneously with the continuity equation and the boundary condition 
equations for the Fourier-Chebyshev coefficients. A nearly block-tridiagonal matrix 
results in the channel problem, in which Cartesian coordinates are used. It was 
found that in cylindrical coordinates a much more complicated matrix results. 

In mother approach, Orszag & Kells (1980) used a fractional-step scheme, using 
Chebyshev polynomials, which seems to be quite efficient for the channel problem. 
Similar schemes have been used in cylindrical coordinates for the flow in a pipe 
(Patera & Orszag 1981) and for Taylor-Couette flow (Marcus, Orszag, & Patera 
1982); however, they result in matrices that are solved in 0(N 2 ) operations, where 
N is the number of Chebyshev polynomials. In the fractional- step scheme used by 
these authors, each time-step is split into three independent “corrections.” First, 
the nonlinear terms are explicitly time- advanced, yielding an intermediate field v 
Then the pressure correction is applied, enforcing the continuity constraint on the 
second intermediate field v. Finally, the viscous correction is performed, allowing 
the imposition of the boundary conditions on the velocity field at the new time-step. 
Note that imposing the continuity constraint on the intermediate field v leads tv. an 
error in the continuity equation of order At/Re for the final field. This appears to 
cause no serious problems in the channel calculations of Orszag and Kells; however, 
Marcus, Orszag, & Patera (1982) experienced some accuracy and stability problems 
related to the splitting when calculating Taylor-Couette flow. 

A third method is given by Kleiser & Schumann (1980), who developed a method 
for the channel problem using Chebyshev polynomials. It is quite similar to the 
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fractional-step scheme mentioned above. However, a procedure is used to deter- 
mine the proper boundary conditions for the pressure to ensuio satisfaction of the 
continuity equation. This method is quite efficient for the channel problem, and 
Marcus (1083) has extended this method to cylindrical coordinates to solve Taylor- 
Couette flow. 

Filially Leonard &; Wray (1982) recently developed a new method and applied it 
to flow in a pipe. A spectral representation based on Jacobi polynomials was used 
which inherently satisfies the continuity and boundary constraints. 

In the method used here, we follow Leonard & Wray (1982) and represent the 
velocity field using vector functions that satisfy the continuity equation and bound- 
ary conditions. In this way, these constraints are automatically satisfied. Satisfying 
the continuity equation also remove® n degree of freedom, and, since the pressure 
is eliminated from the equations, oniy two dependent variables are left. Applica- 
tion of this method in both Cartesian and cylindrical coordinates using Chebyshev 
polynomials is described in §2, and some implementation details are presented in §3. 

1.4 Motivation and Objectives 

The available turbulence data for curved wall-bounded flows is deficient in sev- 
eral ways. First, most of the statistical correlations appearing in the Reynolds-stress 
equations, especially those involving pressure and spatial derivatives, are very dif- 
ficult if not impossible to measure, and thus are unavailable. Second, detailed 
statistical and time-dependent data needed to study the structure of turbulence are 
scarce. Finally, in many experiments Taylor- G or tier vortices have been observed on 
concave curved walls. However, detailed data on the variation of relevant turbulence 
parameters with position in the vortices and on the contribution of the vortices to 
turbulent stresses and higher order statistics are not available. 



Turbulent flow simulation is uniquely suited to provide the information outlined 
above. Moreover, such a simulation provides data relevant to general wall-bounded 
flows since it is anticipated that many of the features of wall-bounded flows are 
qualitatively unaffected by mild curvature. In particular, the structure of the flow 
very near the walls can be studied to discern the curvature-induced differences 
between the concave and convex walls in a curved channel, and, in so doing, the 
near wall features of noncurved wall bounded flows are observed as the features 
common to both walls. A direct numerical simulation of a wall-bounded flow is 
particularly valuable for verifying the results of large-eddy simulations, because 
a direct simulation has not previously been performed of a fully developed, wall- 
bounded flow. 

The study reported here was undertaken to perform such a simulation, and had 
the following specific objectives: 

(i) Develop and verify a fully spectral numerical method for incompressible flow 
between parallel planes and concentric cylinders. 

(it) Demonstrate the feasibility of performing direct simulation of a fully developed 
wall-bounded turbulent flow. 

(tit) Develop a detailed data base for curved turbulent channel flow to be used to 
study curvature effects and the more general features of wall-bounded flows, 
(tv) Assist turbulence modelers by computing quantities of interest that are difficult 
or impossible to measure. 

(v) Investigate the presence and role of turbulent Taylor Gortler vortices in this 
flow. 

(vt) Investigate the structure of turbulence present in curved-channel flow and more 
general wall-bounded flows, using statistical methods and computer visualiza- 
tions. 
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2. Mathematical and Numerical Considerations 


In this section, a method is described for numerically solving the incompressible 
Navier-Stokcs equations in a bounded domain. Much of the material presented 
here appears in Moser, Moin, k Leonard (198.3) and is included in the interest of 
completeness. In the time- advancement procedure used here the viscous term is 
treated implicitly, whereas an explicit scheme is used for the nonlinear (convective) 
terms. In this mixed explicit-implicit time-diffcrencing, the explicitly treated terms 
act as forcing terms to the implicit part of the calculation. In essence, then, an 
implicit time- advancement procedure is needed for the forced Stokes equations (the 
nonlinear term is replaced by a forcing term). It will be convenient in much of 
the discussion that follows to consider only the Stokes equations; however, the 
Navier- Stokes equations can be easily solved with any scheme for solving the Stokes 
equations, given a technique for computing the nonlinear terms, vxu. 

We consider the forced Stokes equations, 

= —VP - ^-VxVxv + f , 
at Re 

V • v = 0, (2.0.1) 

v = 0 at the walls , 

where v is the velocity vector and f is some forcing function. 


2.1 Divergence- Free Vector Expansions 

We seek a finite spatial representation of the velocity vector, v. Since v is con- 
strained to satisfy the continuity equation. and the boundary conditions, we choose 
a representation v«, which inherently satisfies these constraints, 

J 

v 4 (x ,y,z,t) = '52 a i( t ) w A x >v > z ) » (2-i.l) 

j = o 
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where otj(t) arc the coefficients of the expansion, ai:d w j(x,y,z) are a set of vector 
functions chosen to satisfy 

V • w j — 0, w j ~ 0 «*vt the walls . (2.1.2) 

The representation must also be complete so that for sufficiently large J a'l vector 
fields of interest can be represented by (2.1.1). 

The representation (2.1.1) is substituted into the Stokes equations (2.0.1) and a 
weighted residual method is used to obtain ordinary differential equations for the 
coefficients c*j(t). This involves dot multiplying the equations by a set of weight 
vectors and integrating over the computational domain. Vector weight functions 
$;-/ are chosen such that 


V • — 0, • n = 0 at the walls , 


(2.1.3) 


where n is a unit vector normal to the wall. When the weight vectors are formed 
in this way it can be shown, using integration by parts, that the pressure term is 
eliminated from the resulting equations. The result is 


£ TT L • w ’ * = w > 

+ / 

Jd 


dv 


f dv . 


These equations can be written in the compact form, 

where A and B are (J + 1) x (J + 1) matrices with elements 

Aj'j = / fcj'-Wj du, 

Jd 

Bj>,j — ~ I ®j' • Vx Vxw j dv , 
Jd 


(2.1.4) 


(2.1.5) 


( 2 . 1 . 6 ) 
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where a is the.vector with elements ay and F is the vector with elements 


Fy- = / $y. f dv . (2.1.7) 

Jd 

Equation (2.1.5) is a system of linear ordinary differential equations which can 
be solved numerically using any standard time- discretization scheme. It should 
be noted that even an explicit scheme will require “inversion” of rhe matrix A\ 
therefore, unless A is much more sparse than 8, there is no computational advantage 
in using an explicit scheme. 

It will be useful to consider the scheme described above from two ;; mints 
of view. First, it is shown in Appendix A that if (2.0.1) has a solution, it may be 
obtained by solving the following “weak” problem, which is due to Leray (1934): 

v G M = {u G Hj, V • u = 0} , 

J t M = -i(Vxu,Vxv) + (u,f) , (2.1.8) 

VugV, 

where (u,v) is an inner product defined as 

(u,v) = [ u vdv, (2.1.9) 

Jd 

and V is the space of all divergence-free vector functions satisfying the no-slip 
boundary conditions with at least one square integrable derivative . It is clear that 
in the numerical method described above, the functions Wy are a basis for a finite 
dimensional subspace of V (wy has at least one square- integrable derivative), and 
the functions $y< also span a subspace of V. If the numerical method is to be 
consistent with (2.1.8) both the finite dimensional subspaces spanned by wy and 
$y< must converge to “V as J becomes infinite. 

The second point of view is related to the fact that any vector field u can be 
uniquely decomposed into two fields S and V</>, where S satisfies the no-flow-through 



condition and is divergence free (Ch vin & Marsden 1979). There is a projection 
operator P which maps functions u into divergence-free functions S. By rewriting 
(2.0.1) as 

d\ 1 

— + Vi°- - — VxVxv + f (2 1.10) 

at Re 

it can be seen that since d\/dt is divergence-free it is the projection of the right- 
hand side of (2.1.10), that is, 

^(-i-vxvx. + r). 

As is shown in Appendix A, the numerical method described here can be viewed as 
an approximation to this projection operator, P J . Let A be the difference between 
the right-hand side of (2.1.10) and its approximate projection, 

(4 

If the projection were exact, A would be the gradient of the pressure P\ however, 
since the projection is approximate (because of the finite spatial resolution), A is 
not necessarily the gradient of a scalar. The pressure can be defined through the 
exact projection of A, 

VP = A - P{ A) , (2.1.13) 

and P( A) is a truncation error associated with the projection; it will be called the 
projection error. If the numerical method is to converge, the projection error must 
go to zero as J becomes infinite; also, we expect that for problems with adequate 
resolution the projection error will be small. Note that the presence of projection 
error is a manifestation of poor spatial resolution, and it can be monitored to 
ascertain the adequacy of the reso'ution. 


A = -~V xVxv + f — P J 
Re 


VxVxv+f 


)■ 


( 2 . 1 . 12 ) 
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2.2 Projection Error 


The projection operator P discussed above maps any vector function into a vector 
function satisfying the no-flow-through condition, not the no-slip condition. In fact, 
for a general vector field, there is no decomposition into S and V<t> such that S is 
divergence free and satisfies the no-slip condition (Cliorin & Marsdcn 1979). Thus, 
in general, the right-hand side of (2.1.11) is not guaranteed to vanish at the walls. 

Consider the Navier-Stokes equations (f is replaced by v x u in equation 2.0.1). 

1 

— + VP = VxVxv + v x u , (2.2.1) 


or 


dv 

dt 


>(-± 

{ Re 


Vx Vxv + v 


x . 


( 2 . 2 . 2 ) 


If v is such that the right-hand side of 2.2.2 is zero at the walls, v will be said to 

♦ 

be a compatible velocity field. This is the compatability condition pointed out by 
Mom, Reynolds, & Ferziger (1978). It has been shown (Temam 1983) that with 
an initial condition that has n > 2 square-integrable derivatives, the Navier-Stokes 
equations have a unique solution for 0 < t < T for some finite time T, which has 
n square integrable derivatives and is compatible in the above sense. Also, dv/dt 
has, in general, n — 2 square-integrable derivatives. 

At best then, dv/dt has two fewer square-integrable derivatives than v; at worst 
v may not be compatible so that dv/dt does not have even one square integrable 
derivative. In either case, dv/dt is not as smooth as v and is, therefore, more 
difficult to represent with a numerical method. The spatial discretization error 
in dv/dt may thus be much larger than in v, and this error contributes to the 
projection error in the method presented in §2.1. 

The computations of the various states of Taylor-Couette flow described in §4 
have been examined for projection error. In the low Reynolds-number cases (Taylor 
vortices and wavy Taylor vortices), where the velocity field and dv/dt are very 






smooth, the relative projection error |P(A) 2 |/|A 2 | (here j • j signifies the average 
over a surface parallel to the walls) was less than 10~ 6 everywhere. In the higher 
Reynolds- number case of modulated wavy vortex flow, the velocity field and par- 
ticularly the time-derivative field are less smooth and there is some evidence of 
weak turbulence. For this case the the relative projection error was of the order of 
10~ 5 , except very near the wall where it was 10~ 2 . These results confirm that the 
projection error is related to the lack of smoothness of dv/dt. 

Projection errors have also been observed in the turbulence calculations being 
reported here. Values of |P(A) 2 |/|A 2 | as large as 1 have been observed near the 
walls. This relative error falls off quickly to less than 10 -2 going away from the 
wall. Values of ||F(A)|| 2 /||A|| 2 (|| • || signifies the L 2 norm) and ||^(A)||/jjc , v/3t|| 
were of the order of 10 ~ 3 . It is not surprising that the error is so much larger 
at the walls than away from the walls in view of the discussion above. To inves- 
tigate these errors further, the approximate projection P 2J (using twice as many 
expansion functions) was performed and compared with P J reported above. The 
projection error P( A 2J ) at the walls was reduced by an order of magnitude from 
the value reported above, and the errors in the center of the channel were reduced 
by two orders of magnitude. The reduction of error at the wall when the more 
accurate projection is used indicates that most of the error in P J at the wall is not 
caused by the incompatibility of the velocity field. If the velocity field were truly 
incompatible the wall error would not be reduced with improved approximations to 
the projection. Instead, the large values of relative projection error at the walls is 
caused by the finite spatial resohrtion used in the projection. The more rapid re- 
duction of projection error away from the wall is further indication of the difficulty 
of representing dv/dt near the wall. These projection errors are indicative of the 
marginal resolution of the computations. 
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2.3 Plane Channel 


In order to use the method described in §2.1 the expansion functions and weight 
functions must be chosen. In this section, functions for plane-channel flow are 
presented. Although plane-channel flow is not under investigation in this work, it 
is instructive to Grst present the numerical method for this simple geometry. The 
functions used in the curved channel will be presented in §2.4. 

The plane-channel flow is homogeneous in the two directions parallel to the walls 
(i and z). Therefore, in these directions, periodic boundary conditions and Fourier 
expansions are used. Only the expansions in the direction normal to the walls ( y ) 
need to be determined. It is assumed that the walls are located at y = ±1. The 
expansion functions of (2.1.1) become 

w jmi = Uj(y,k x ,k z )e' k * x e tk ‘*, (2.3.1) 

where 

/9t rnO^r 

k x = — , — L<1 < L\ k z = - M < m < M (2.3.2) 

L x L z 

are the wave numbers; L x , L z are the periods in x and z , respectively; L, M are half 
the number of Fourier modes used, and Uj(y;k x ,k z ) are a set of vector functions 
chosen such that 

V • w j m i = 0 and u,(y = ±l) = 0. (2.3.3) 

The weight functions are chosen as 

= *Ay\k x ,k z )e~ ik '* (2.3.4) 

with 

V • = 0, and ®y*(y = ±1) • n = 0 , (2.3.5) 

where n is a unit vector normal to the wall. After using the orthogonality property 
of the complex exponentials to evaluate the x and z portions of the integrals in 
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(2.1.4), the weighted residual method of §2.1 yields the following s< t of equations 
for each wave-number set k x ,k x : 



(2.3.6) 


where .and Uy depend parametrically on k x and k M , f is the Fourier transform 
of f, and V x is the Fourier transformed curl operator. Note that the subscripts l 
and m were introduced for the expansions in the x and z directions but have now 
been dropped for brevity. 

We now turn to the problem of choosing the vectors 'Py* and Uy . There is consid- 
erable freedom in this choice. The vectors presented here were constructed to yield 
matrices A and B in (2.1.5) which are banded with small bandwidths. This is easily 
accomplished when one of the wave numbers (say k z ) is zero. First, vectors will be 
constructed for this special case; later the results will be extended to the general 
case. 

It becomes necessary to split the sets of vectors and u into two classes, (\P + , \& _ ) 
and (u + ,u~), respectively, each class having a different functional form. This is 
equivalent to independently representing two components of the velocity vector, 
with the third determined by the continuity equation. To obtain tightly banded 
matrices, we desire that the equations for u + be decoupled from the equations for 
•i~ , that is, 


r 1 

r l 

/ 9+ • u: dy = 0, 

/ 9 1 , • V x V x u ~ dy = 0, 


J - 1 

r 1 

pi 

/ 9 7, . ut dy = 0 , 

1 • Vx Vxuy dy = 0 

J-i 1 

J-i 
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A convenient choice which satisfies this requirement (for k t = 0) is 



(2.3.8) 


where g, Q, h, and P are indexed functions of y to be chosen later, and the 
superscript prunes (not on subscripts) indicate differentiation wkh vespcct to y. 
can be easi'y verified that these vectors satisfy the continuity conditions (2.3.3) and 
(2.3.5); the boundary conditions will be satisfied if 


0j( y = ±l)=0, hj(y = ±1) — 0 , 

g'A y = ±l) = 0, Qj.( y = ±l)=0. 


(2.3.9) 


It will also be useful in reducing the band widths of matrices A and B , to require 
that 


Pp(y = ±l) = 0 and Qj’{y = ±1) = 0 . (*.3.10) 


The integrals in (2.1.6) can be evaluated using the identity in Cartesian coordi- 
nates (for V • u = 0), 


Vx Vxu = — V 2 u, 


(2.3.11) 


and integrating by parts with the results 

Ah = /_* n ■ *■; <hi = - /‘ QAh.) dy . 

Bpj = - f *>. -^xvxu; <(» = -/' (CQAiCii) iy , 

/’ (2.3.12) 

A?J = J _ x *;■ • »,• dy = Jf_ i P,.(£M iy . 

SJ.J = -J VxuJ iy = / Pj'(Chj) dy , 
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where £ is the Fourier- transformed Laplacinu operator, 

£ - — - fc 2 
^ “ dy 2 x ‘ 


(2.3.13) 


Therefore, when the vectors (2.3.8) are used, (2.3.6) is decomposed into two inde- 
pendent equations for the coefficients of u f and u~, 




(2.3.14a) 

(2.3.14b) 


with matrices as defined in (2.3.12) and F + and F" defined a3 follows: 

Fp = j l {k x Q } ,i y -iQpi x )dy, 

F ? = j x p j'f* d y • 


(2.3.15) 


It should be mentioned that in Cartesian coordinates (2.3.14) can be derived in 
a more straightforward way (this will not be the case in cylindrical coordinates; see 
§2.4). Equation (2.3. 14b) is readily obtained from the z equation of (2.0.1), after 
Fourier transforming (again fc, =0), 


If the representation 


di - 1 0* -l f 

3T~Re 


v*=Yi a i h ’ 


(2.3.16) 


(2.3.17) 


is used in a weighted residual method with weights Py, (2.3.14b) is obtained. For 
(2.3.14a), the curl operator is used twice on (2.0.1), which yields 


'l(VxVxv) = -J-(Vx) 4 v + VxVxf . 
at Re 


(2.3.18) 
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After Fourier transforming and using (2.3.11), the y equation of (2.3.18) is the 
fourth order equation 

= -~U\) + x Vxf)„ . (2.3.19) 

The representation 

= (2.3.20) 

j 

and the weight functions Qj> yield (2.3.14a). The continuity equation evaluated at 
the wall requires 

A 

^(y = ±l) = 0, (2.3.21) 

which provides the additional boundary conditions on v v to make (2.3.18) well 
posed. The x velocity is determined from the continuity equation. Thus, for this 
case, solution of (2.0.1) is equivalent to the solution of (2.3.16) and (2.3.19). These 
equations were solved in a method used by Patera & Orszag (1981). 

Extension of the vectors used above to the general case when k x ^ 0 and k z ^ 0 
is straightforward. By rotating the coordinate system about the y axis, the general 
problem can be reduced to the k 2 = 0 case already discussed. The new axes (x 1 
and z') are rotated such that the x' axis is aligned with the vector 

k x e x + k z e z , (2.3.22) 

where e x and e* are unit vectors in the x and z directions. Then, 

k x > = (kl + k\yi 2 , k t , = 0, (2.3.23) 

and the vectors (2.3.8) can be used. This is the coordinate transformation at the 
heart of Squires’ theorem in the hydrodynamic stability of plane parallel shear flows 
(Stuart 1963). 
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Finally, the .vectors defined in (2.3.8) are incomplete when both k x and k z are 
zero. For this case, the following set of vectors is used: 


u y “ 

fhj) 

0 

ii 

1 

3 

o o 


\°J 

f p r > 


° \ 


0 

* = 

0 


^ 0 ; 


{ p j’j 


(2.3.24) 


This leads to identical matrices for the plus and minus equations; the derivation 
follows that for (2.3.14b). 

The expansion functions g, Q, h, and P must now be chosen. Strictly orthogonal 
functions (which would lead to diagonal matrices A and B in 2.3.14) should not be 
used because requiring orthogonal functions to satisfy boundary conditions (2.3.9) 
imposes extraneous conditions on higher derivatives of the functions, which degrades 
the rapid convergence of the method (Gottlieb Sc Orszag 1977). Instead, we use 
quasi- orthogonal functions, which lead to banded matrices A and B, and do not 
suffer from this convergence problem. Quasi-orthogonal functions are constructed 
from a set of orthogonal functions which admit general boundary conditions (see 
Gottlieb Sc Orszag 1977 for a discussion of this class of functions). Since these 
functions do not inherently satisfy any boundary conditions, boundary conditions 
are imposed by forming linear combinations of them to make the quasi-orthogonal 
function. This construction must be done in such a way as to make matrices A and 
B , which involve integrals of the functions and their derivatives, banded. Orthogonal 
polynomials are suitable for this purpose because they satisfy rectirsion relations 
which make this construction easier. 

The Chebyshev polynomials have been chosen for this application because they 
have two properties that are particularly useful, (t) They are related to the cosine 
function through a coordinate transformation (Fox Sc Parker 1968); this allows the 
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use of the fast Fourier transform algorithm in evaluating F. (u) They arc partic- 
ularly efficient for resolving boundary layers near the walls (y ~ ±1) (Gottlieb Sc 
Orszag 1977). The functions g, Q, h, and P are constructed so that they and their 
derivatives have simple forms when expressed as linear combinations of Chebyshev 
polynomials, which guarantees that the matrices A and B are banded. This can be 
accomplished by manipulating the recursion relations for the Chebyshev polynomi- 
als. The constructed functions g, Q, h, and P are, 


9 i = (i - y 2 ) 2 ?V(y), 


>-(£ 


T i+ 2 


2T, 


+ 1 ) 0 + 1)0 

P i = (T,_ 1 -T y+1 )/2/( l-y 2 ) 1/2 


hj = (1 - y 2 )Tj(y), 

h) + M)/ 4(1 - s3)1/2 ’ (2 - 325) 


where Tj is the Chebyshev polynomial of order j and the factor 1/(1 — y 2 ) 1 ^ 2 ap- 
pearing in Q and P is the weight for the Chebyshev orthogonality relation. 

Other orthogonal polynomials can be used to construct quasi-orthogonal func- 
tions. This is a consequence of the recursion and differential relationships that 
orthogonal polynomials satisfy (the Chebyshev relationships are particularly sim- 
ple). Thus, there are many possible sets of quasi-orthogonal functions that can be 
used to meet requirements that might be imposed in a given problem. For example, 
Legendre polynomials could be used in this problem instead of Chebyshev polyno- 
mials; in that, case the g and Q functions would be identical, as would the h and P 
functions. The resulting approximation would have identical expansion and weight 
functions so that the analysis in Appendix A would apply. A second example is the 
functions based on shifted Jacobi polynomials used by Leonard Sc Wray (1982) in 
the calculation of flow in a pipe. 

With the choice of the functions g , h, Q and P the method is completely defined. 
Of great interest is the amount of computation required to implement this method. 
After taking advantage of the decoupling f even and odd functions, both matrices 
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A' and S r are banded with two subdiagonals and four superdiagonals, and A" and 
B~ have one subdiagonal and two superdiagonals. Also, the elements in the matrices 
depend only on the square of the wave numbers k x and k t so that at least two wave- 
number sets can be solved simultaneously. The result is that for each wave-number 
set k x , k z , the matrices arising from the implicit time- differencing of (2.3.14) can 
be solved in 30 N real additions and multiplications, where N is the number of u f 
vectors used (three less than the highest-order Chebyshev polynomials used). There 
is some additional computation required to calculate the forcing vector F and to 
perform the coordinate rotations discussed earlier. The total cost is then about 
50/V operations. This compares favorably with the 376 N operations required for 
matrix solution in the finite difference method used by Moin & Kim (1982). Thus, 
the method is operationally efficient, in addition to offering savings in computer 
storage by reducing the number of independent variables. 


2.4 Curved Channel 

In this section, expansion and weight functions to be used for a curved charnel 
will be described. We shall consider the flow in an annulus. In this case the flow is 
assumed to be periodic in the axial ( z ) and the azimuthal (0) directions. The inner 
and outer walls are located at r = r, and r = r 0 , respectively. Using representation 
and weight vectors as in (2.3.1) and (2.3.4), a result (for cylindrical coordinates) 
similar to (2.3.6) is obtained. 

t % r *’■ ■ % * * = -b t jp ■ ♦* w * 

i-° r% >=0 r ‘ (2.4.1) 

r° 

+ / V]' ■ f r dr , 

Jfi 


where ®yr(r) and u, (r) depend parametrically on k$ and k z , which are the 6 and z 
wave numbers. 
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Again, it is .desirable to choose $ and u to make A and B tightly banded. The 
first step is to attempt to uncouple the equations for the two sets of vectors 3» 1 ', 
u + , and \P“, u“. In cylindrical coordinates, however, the appearance of uo in the 
r-component of V x V x u (for V • u = 0) and u r in the 0-component makes the 
decoupling more difficult than in the Cartesian case. 

The following vectors satisfy the decoupling requirement, though they have an 
important defect to be discussed later. 



= V*xV*x 


u; =Vx 


$: ( =V*xV*x 



f ~ ik x Qj 
- ki9j 
9- + ^ 


9j 


f-w) 

f - ik z Qj \ 

9i = 

k *9j 

l 0 J 




(2.4.2a) 


(2.4.2b) 


(2.4.2c) 


(2.4.2d) 


where Qj> and py are indexed fimetions of r (not the same as in (2.3.25)) and 
V* x is t^e complex conjugate Fourier-transformed curl operator. Satisfaction of 
the continuity equation is guaranteed by the identity V • (V x u) = 0. In order to 
enforce the boundary conditions, we require 

9j{r ~ r *> r °) = 0 . Qj'{r - r i,r Q ) = 0 , 

9j{r = r,-,r 0 ) = 0 , <?'. (r = r„ r 0 ) = 0 . 

Equation (2.4.1) now decomposes into two equations, as in (2.3.14). However, the 
vector functions defined in (2.4.2) are an incomplete set, because the imposition of 
the boundary conditions in (2.4.3) forces the 0-velocity to satisfy the condition 

du 0 


(2.4.3) 


dr 


= 0, at r = r,-,r 0 , 


(2.4.4) 
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f 0 \ 

( 0 ’i 

u» = 

- 2k z hj , 

®J. = - kyp,. 


^ (2 k 0 /r)hj J 

^ kgrPj / ; 


which is too restrictive. To alleviate this problem, the vectors in (2.4.2) are aug- 
mented with two extra vectors, 


(2.4.5) 


where j and j' can be 0 or 1. With this addition, dug/dr can have arbitrary 
values at the walls. It can be shown that a finite number of expansion functions 
u + , u~, u° can exactly represent the Fourier components of any divergence-free 
velocity field with polynomial dependence on r that satisfies the no slip boundary 
conditions; therefore, the expansion functions are complete. When the extra vectors 
arc included, equations (2.3.14) become 

a: 


dt Re’ 

" da = + S 0 a") + F-, 


dt 


A *w +A °-^r +A ’^ = + flv + s >") + F " 

where the matrices are defined atf 


(2.4.6a) 

(2.4.6b) 

(2.4.6c) 


W)j'j = f • u? r dr, 

Jri 

I'" • VxVxuy r dr. 

Jti 


(2.4.7) 


The subscripts/superscripts 7 and /? can be the symbols +, -, or 0 as used above. 
The equations in (2.4.6) axe coupled, but in : weak way that will not significantly 
affect the computational eflciency of the method (see below). 

The functions g, h, Q, and P are again constructed with Chebyshev polynomials, 

Qj = r ( 1 - v 2 ) 2 Tj{v) » h i = K 1 - y 2 )^(y) > 

1 j-f 2 2T, 


(ifi 


+1) (j+i)(y 

p i = (Tj-i - T j+I )/2l(l - y 2 ) 1/! , 


h) + & (2 - 4 - 8) 
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where y = (2 r -r D - ri)/(r 0 - r,), so that y is --1 where r =- r, and +1 where 
r = r„. These are the same functions used in the Cartesian problem, except for the 
factors of r and These factors are included to cancel the various 1/r’s appearing 
in the C operators, which is necessary if the Chebyshcv orthogonality relations are 
to be used to evaluate tie integrals in (2.4.7). 

The coupled equations (2.4.6a-c) can be written as a single equation, as in (2.1.5) 
(a is composed of the vectors or + ,ar, and a°). The resulting matrices A and B 
have the special form shown in Figure 2.1. Also shown are the bandwidths for 
the various nonzero bands in the submatrices. Though this matrix is not strictly 
banded, it can be solved with no difficulty. As in the Cartesian case, there are wave- 
number synur etries that allow the solution of more than one wave-number set at a 
time. Including these symmetries, the operation count for the matrix solution for 
each wave- number set k$, k r is 235 AT additions and multiplications, where N is the 
number of u + vectors (four fewer than the highest order Chebyshev polynomials 
used). 

The representation presented above is incomplete when k g = 0. The following 
vectors are used in the special case k x = 0, k u 0: 



where j, h, Q , and P are as defined in (2.4.8). When k g = 0 and = 0, the 
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vectors in (2.4.9) are incomplete'. 


In this case the following vectors arc used: 


u 





(2.4.10) 


The vectors in (2.4.9) and (2.4.10) are very similar to the ones used in the Cartesian 
case and lead to uncoupled sets f equations. Solution for these cases requires less 
computation than the general case considered above. 
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3. Implementation Details 


The numerical method outlined in §2.4 has been implemented on the Cray-XMP 
computer. The code CURVE which implements the method was written in VEC- 
TORAL, which is a language developed at NASA- Ames Research Center for large- 
scale computations on vector processors (Wray 1983). Some of the important fea- 
tures of CURVE and some of the other codes used in conjunction with it are dis- 
cussed in the following sections. 


3.1 Time Advancement 


The method outlined in §2.4 is used for the spatial discretization of the Navier- 
Stokes equations. For the tiine-discrctization, the viscous terms are advanced us- 
ing the Crank-Nicholson scheme (trapezoidal rule), whereas second order Adams- 
Bashforth is used for the nonlinear terms (f = v x u). The resulting equations 
are 


{* - ^ B ) “" +l - {* + £,») + f < 3F ” - F "~‘> ■ < 3U ) 


where 



tyj, • (v n x w n ) rdr, 


(31-2) 


A, B, and a are a.- defined in §2.4 and the superscripts denote time- level. Each 
time-step consists of five steps: 

(») Compute v x u from v n and u n 

(it) Compute F n , save it, and add ~*F n to (A + - ^F n-1 saved from 

the previous time-step 
(tit) Solve 3.1.1 for a n+1 

(tv) Compute (A + ^B)a n+1 and add it to -^F n saved from step (it) 

(v) Compute v n+1 and w n+1 from a n+1 
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Iu the first ^tep, v x u is computed by evaluating v and w on a grid of quadra- 
ture (collocation) points and performing the required multiplications. To eliminate 
aliasing errors, the number of quadrature points in each direction is chosen to be 
3/2 the number of modes in that direction (Paterson & Orszag 1971). To evalu- 
ate v and u on the grid and to express v x w in terms of Fourier functions and 
Chebyshev polynomials the fast Pourier-transform algorithm is used. Once an ap- 
proximation for v x w in terms of Fourier functions and Chebyshev polynomials is 
obtained, the integrals in (3.1.2) are performed semi-analytically to obtain F. This 
is accomplished by using the analytically prc-computed values of the integrals 

f ° Tj* y ■ e< rdr, (3.1.3) 

Jr, 

where e t - are the unit vectors in each of the coordinate directions. Thus, the evalu- 
ation of the integrals in (3.1.2) reduces to a matrix multiply. This matrix is banded 
for 9j, as defined in (2.4.2), (2.4.9), and (2.4.10). 

Solving (3.1.1) for a n+1 once the right-hand side is computed requires the solu- 
tion of a system of linear equations defined by the matrix [A — ^^B), which has the 
banded structure shown in Figure (2.1). This is accomplished using Gauss elimina- 
tion without pivoting. The elements of A and B axe pre-computcd by analytically 
evaluate .ig the integrals in (2.4.7). Once a n+1 is known, computing (A + 2 ^B)<* n+1 
is a handed matrix multiply. 

To compute v x u we must obtain the velocity field v expressed in terms of 
Chebyshev polynomials instead of the divergence free vector representation. Since 
the vector basis functions u* defined in (2.4.2), (2.4.9), and (2.4.10) are known 
in terms of Chebyshev polynomials, this is straightforward and again involves a 
handed matrix multiply. Once v is known in terms of Chebyshev polynomials and 
Fourier functions, u) is obtained by differentiation. 

The various matrices required in the above computations axe computed once and 
stored to be used at each time-step. However, the matrices are functions of the wave 


31 



numbers ko and k x . To avoid storing all these trices for each wave number, they 
are decomposed into sums of several matrices with powers of the wave numbers as 
coefficients; for example 


dJ 


A(ko, k x ) = Ai + k'fiAi + k 2 x Az . 


(3.1.4) 


where Ai, A?, and ^3 arc constant matrices, independent of kg and k x . The com- 
ponent matrices are stored, and A(kg,k z ) and the other required matrices are re- 
constructed as needed. The component matrices are evaluated using a FORTRAN 
routine (GENMAT) which evaluates the class of integrals 



° d ' [ r *Q>'( r )] 

dr' 


dr 1 


r n dr , 


(3.1.5) 


analytically given the representation of gj and its derivatives up to order l in terms 
of Chebyshev polynomials, and the representation of Qj> and its derivatives up to 
order t in terms of Chebyshev polynomials divided by (1 — y 2 ) 1 / 2 . GENMAT uses 
the Chebyshev orthogonality relation to perform these integrals. 


3.2 Data Management 

The turbulence simulations being reported here were performed with 65x 128x 128 
(~ 10 6 ) Fourier/Chebyshev modes. Because the amount of data involved is much 
larger than will fit in the central memory of the Cray-XMP, it is necessary to 
maintain the data on a secondary storage device and perform the computations 
by loading small sections of the data at a time into central memory. To avoid the 
relatively long access times inherent in disk-based mass storage, a very large (16 
million 64-bit words) solid-state memory called the SSD (Solid State Device) is used 
for the secondary storage. In addition, to reduce the amount of storage required 
on the SSD, the data are packed from 64-bit words to 32-bit words, thus reducing 
the storage and precision by a factor of 2. Note that the 32- bit precision is used 
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only on the SSD; all computations are done with 64- bit precision. This reduction in 
precision each time data is stored on the SSD has minima! impact on the accuracy 
of the computations. 

To use secondary storage effectively, the computations and the data base must 
be carefully structured. For each time step, two passes are made through the data 
base (see the flow chart in Figure 3.1). In the first pass (PASSl), vxuis computed 
and the required Fourier transforms in the 0 and z directions are performed. In the 
second pass (PASS2), the Chebyshev transforms in the r direction are performed 
as arc the matrix operations outlined in steps (it) through (tv) in §3.1. On leaving 
PASS2 and entering PASSl, the main data base (called VDATA) must contain the 
three components of velocity, as well as dv 0 /dr and dv z /dr. These two derivatives 
are needed to compute the vorticity in PASSl, but they must be computed in 
PASS2. On leaving PASSl and entering PASS2, the VDATA data base contains 
the three components of v x u; the space that previously contained the derivatives 
is empty. The data in VDATA are Fourier transformed in the 9 and z directions, 
but are not Chebyshev transformed in the r direction; for example, the velocity is 
stored in VDATA as v(r, ko,k z ). In addition to the data in VDATA, the result of 
step (tv) in §3.1 is computed in PASS2 and must be saved to be used in PASS2 
on the next time-step. These data are stored in the secondary data base, RDATA, 
which is accessed only in PASS2. 

When performing the Fourier transforms, all the data to be transformed must be 
in core memory. Therefore, in PASSl, where the transforms in 9 and z are done, 9-z 
planes of data at individual r locations must be brought into core. In PASS2, we 
only need to have all the data for a given wave-number pair kg, k z in core; however, 
data were loaded into core in r-0 planes to facilitate vectorization. To allow the 
data in VDATA to be brought into core in either r-9 or 9-z planes it is divided 
into an array of “drawers,” each of which can be accessed separately. Each drawer 
contains data for all values of kg and for a few points in r and a few values of k z 
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(5 and 4, respectively, in the current calculations). To bring in a 0-z plane, a row of 
drawers at constant r is accessed; to bring in an r-0 plane, a column of drawers at 
constant k z is accessed. Each drawer is arranged to contain the data for both wave 
numbers k z and -k z . This is advantageous because in PASS2 the matrices to be 
solved for k z and — k z are identical. The data contained in the RDATA data base 
are only required in PASS2, so RDATA need not be structured in drawers. It is 
divided into slabs that contain several r-9 planes of data. Each slab contains data 
for the same values of k z as a column of drawers in the VDATA data base. See 
Appendix B for more information on the data bases. 


3.3 Statistics 


A variety of statistical correlations have been collected from the turbulence com- 
putations. Among the statistics gathered were mean velocity, turbulent intensities, 
rms vorticity fluctuations, velocity skewness and flatness factors, and Reynolds shear 
stress. All of these quantities were calculated as a function of r by averaging the rel- 
evant quantities in the 0 and z directions and time. In addition, the mean velocity 
and Reynolds-stress tensor were obtained as a function of r and z by averaging only 
in 9 and time. To improve the statistical sample, these two-dimensional statistics 
were also averaged over the mirror image of the computed flow ( z is mapped to 
—z). This was done because the Navier-Stokes equations are invariant with respect 
to a reflection so that the mirror-image flow is an equally valid solution. 

All terms in the Reynolds-stress tensor balance equations have been computed. 
Again, the relevant quantities were averaged in the 9 and z directions and in time. 
Some of these terms involve the pressure, which is not directly available from the 
computations. To obtain the pressure, the following Poisson problem was solved: 


V 2 P =V • A , 

VP • n =A • n at the walls , 


(3.3.1) 
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where n is the unit vector normal to the wall and A is defined in (2.1.12). The 
solution was obtained using a Chebyshev weighted residual method with expansion 
functions satisfying the derivative boundary condition and weight functions Pj> ns 
defined in (2.4.8). The contribution of the projection error to the Iteynolds-stress 
balances was also computed. 

Two point correlations and spectra were compute 1 in 0, z, and time (t); one- 
dimensional and two-dimensional correlations have been obtained. Correlations 
and spectra not involving time were found at selected r locations by computing 
the two-dimensional (0 and z) energy spectra at those locations. These spectra 
were sampled every 20 time-steps throughout the computations to ensure a good 
statistical sample. The two-point correlation functions were obtained by Fourier 
transforming the spectra; they include the averages over the mirror image-flow. 

Frequency spectra were obtained by temporally Fourier transforming the velocity 
field at selected r locations. This Fourier transform was performed as the calcula- 
tions progressed. The discrete temporal Fourier transform is defined as 

X>(* At )«"'“'. (3-3.2) 

v Jfc=0 

At each time-sample point, the velocities in the chosen plane were multiplied by e ,u ' jt 
for each uij. The result was added to the sum of previously calculated products, 
thus producing the sum in (3.3.2). However, it was not possible to use equally 
spaced u>j, as would normally be done in a discrete Fourier transform, because for 
the range of frequencies of interest, this would require too much storage. Instead, 
the Wj's were chosen in a geometric series, so that log(coy) are equally spaced. This 
has the advantage that more data are obtained at the low frequencies where most 
of the energy resides. To obtain the two-point correlation, the frequency spectra 
are Fourier transformed; this was accomplished by first interpolating the spectra 
available at the chosen u>, points to a set of uniformly spaced points w* using a 
cubic spline. 
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4. Verification 


The major verification tests for CURVE have been the computation of several 
states of Taylor-Couette flow. Taylor-Couette flow is the flow between concentric 
cylinders with the inner cylinder rotating. As the Reynolds number increases this 
flow undergoes a series of three transitions to increasingly complicated laminar 
states before the transition to turbulence. These transitions and the intermediate 
states have been under extensive theoretical, experimental, and numerical inves- 
tigation; therefore, there is a large body of information concerning this flow with 
which the numerical results can be compared. 

The first transition is from the one-dimensional laminar Taylor-Couette flow to 
a system of counterrotating, axisymmetric toroidal vortices (Taylor vortices, after 
G. I. Taylor who first predicted their existence, Taylor 1923). In addition to the 
formation of the vortices, this transition is marked by a sudden change in slope of 
the torque required to drive the inner cylinder as a function of Reynolds number. 
Both the critical Reynolds number for this transition and the torque for a Reynolds 
number greater than critical have been computed for two geometries: with inner to 
outer radius ratio rj = ri/r 0 — 0.95 and with 77 = 0.5. These computations were 
performed with no aliasing errors, using nine Fourier modes in the z direction and 
11 Chebyshev polynomials in the r direction. The length in the z direction over 
which periodicity was imposed was chosen to be the wavelength of Taylor vortices 
corresponding to the minimum critical Reynolds number (A/2 6 = 2.009 for the 
narrow-gap case and 1.988 for the wide-gap case; DiPrima & Eagles 1977). The 
results of these calculations are summarized .in table 4.1. 

The critical Reynolds number for transition to axisymmetric Taylor vortices was 
determined by searching for the Reynolds number at which a small disturbance 
neither decayed nor grew. A disturbance of wavelength A was added to the laminar 



TABLE 4.1 


Results of Axisymmetric Calculations 



Critical Reynolds No. 


Torque, G 


Stability 

analysis* 

Present 

calculation 

Re 

Experi- 

mental* 

Present 

calculation 

Narrow gap 
T) = 0.95 

184.99 

185 

195 

5.26 x 10 5 

5.42 x 10 5 

Wide gap 
r) = 0.5 

68.19 

68.2 

78.8 

1.479 x 10 3 

1.487 x 10 3 


t DiPriina & 


Eagles (1977), 


t Donnelly &; Simon (1960) 


solution, and the time- evolution of the first Fourier coefficient was monitored. The 
disturbance would decay rapidly at first, until it consisted of only the least stable 
eigenmode; it would then either slowly grow or slowly decay, depending on whether 
the Reynolds number was above or below critical. Critical Reynolds numbers found 
in this way are presented in table 4.1; note that they are in excellent agreement 
with the analysis of DiPrima & Eagles (1977) for both narrow-gap and wide-gap 
problems. 

For the torque calculations, a disturbance to the laminar solution of wavelength 
A was allowed to grow to steady state. The nondimensional torque G (torque per 
unit length normalized by pi/ 2 , where u is the kinematic viscosity of the fluid) was 
computed from the formula 

G = 27rrRe ^^(rvo) - , (4.1) 

where the overbar denotes average over z. These calculated torques are presented 
in Table 4.1 together with the data of Donnelly & Simon (1960). However, there are 
two reasons why comparisons with the experimental data should be made with some 
caution. First, the axial wavelength of the Taylor vortices in the experiment was 
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not measured. The wavelength corresponding to minimum critical Reynolds number 
was used in the calculations as a good guess, because the Reynolds numbers are 
not far above the critical Reynolds numbers, raid the experimental conditions were 
obtained by slowly increasing the speed of the inner cylinder from zero. Second, the 
experimental torque values in the subcritical regime arc not in very good agreement 
with the torques predicted for laminar Taylor-Couette flow. This is especially true 
for the narrow-gap experiment, in which the experimental torque is consistently 3% 
below the theoretical value. For the wide-gap experiment, the data are within 0.6% 
of the laminar torque for Reynolds numbers far below critical. In light of these 
considerations, the agreement of the present calculations with the experimental 
data of Donnelly and Simon is as good as can be expected (within 3% and 0.5% for 
narrow and wide gaps, respectively). In Figure 4.1, contours of the secondary flow 
stream function are plotted for the narrow-gap case, showing the familiar Taylor 
vortices. 

The next transition is to nonaxisymmetric wavy Taylor vortices. Here, Taylor 
vortices develop waviness in the azimuthal direction, with a range of possible wave 
lengths. These waves travel about the cylinder at about half the speed of the inner 
cylinder. 

The critical Reynolds number for the transition to nonaxisymmetric flow for the 
case A = 2.007, rj = 0.877, and m = 1 (where the ^-wavelength is 27r/m), was 
determined as before by introducing a disturbance and allowing it to grow or decay. 
In this case, the base flow was Taylor vortices calculated with nine Fourier modes 
and 11 Chebyshev polynomials. For the three-dimensional calculations, nine Fourier 
modes were used in the 0 direction. The critical Reynolds number was found in this 
way to be 130, in good agreement with the value of 131 reported by Jones (1981). 
The growth rate and fundamental frequency (precession speed), for an unstable 
nonaxisymmetric mode (Re = 177.6, A = 2.007, r\ = 0.877, m = 6) were found 
by allowing a small disturbance to develop until it consisted of only the unstable 
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mode. Resulting values for growth rate (1.11 X 10~ 2 ) and precession speed (2.54) 
me within O.G% and 0.5%, respectively, of the values reported by Jones (1081). 

Fully developed wavy vortices (Re - 458., A ~ 3.0, r? -- 0.8G8, m — G) were 
computed using 31 Fourier modes in the z direction, 15 in the 0 direction, and 
Chebyshcv polynomials through order 32. Figure 4.2 shows contours of axial veloc- 
ity at r — 0.882r o (close to the outer cylinder), showing the wave in the 0 direction. 
The calculated angular wave speed was 0.3344(1 j, where (7* is the rotation speed of 
the inner cylinder. This is in very good agreement with the experimental value of 
0.334711, of King et a'.. (1983) and the numerical value of 0. 3344(1, obtained with 
a different method using similar resolution (King et al. 1983). 

The final transition studied here is from wavy Taylor vortices to modulated wavy 
vortices. Modulations appear on the azimuthal traveling wave. These modulations 
may be interpreted as a second traveling wave at a different speed. The second wave 
also has a variety of possible azimuthal wavelengths, which give rise to different 
modulation patterns. Modulated wavy vortex flow has been computed for the case 
Re = 1300, A = 2.36, m i = m 2 = 4, rj = 0.877. The precession speed of the 
two waves were computed to be 0.327(1, and 0.437(1,, in good agreement with the 
experimentally determined values 0.33(1, and 0.44(1, (Gorman h Swinncy 1982). 
In addition, it was found that the wave speeds vary by about 7%, depending on 
the position in the modulation cycle. This frequency modulation is also in good 
agreement with the experiments of Gorman and Swinney. Figure 4.3a and 4.36 
show contours of axial velocity at r — 0.886r o (close to the inner cylinder) at times 
separated by 180° in the modulation cycle. Note the more pronounced waviness in 
Figure 4.36 than in Figure 4.3a. 

All the results quoted above show remarkably good agreement with theoretical, 
experimental, and other numerical results. The results of these test cases instill 
confidence in the accuracy of the numerical method and the reliability of the code 
that implements it. 
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5. Curved Turbulent Channel Flow 



In this section, results of the direct numerical simulation of curved, turbulent 
channel flow are presented. These computations were performed using .he numerical 
method described in §2 and the code CURVE described in §3. The curvature 
parameter <5/r c , where 6 is the channel half-width and r c is the radius of curvature 
measured at the centerline, was chosen to be 1/79 = 0.0127. This is within the 
range described by Bradshaw (1973) as mild curvature (6/R ss 0.01); Bradshaw 
suggested that studies on curvature effects should concentrate on mild curvature 
because in problems of aerodynamic interest streamline curvature is most often 
mild. The computational domain (a scale drawing appears in Figure 5.1) has a 
length of 4 tt 0/3 in the spanwise ( z ) direction and subtends an angle of 0.16 radians 
in the streamwise (0) direction, which yields a length of 12.640 along the centerline. 
As will be shown in §5.5 the size of the computational domain was sufficiently large 
for the periodic boundary conditions used in the streamwise and spanwise directions 
to cause minimal distortion of the results. 

Unless otherwise stated, results presented throughout this section will be nondi- 
mensionalized with the shear velocity and channel half-width. However, because a 
curved channel is not symmetric with respect to the channel centerline, the defi- 
nition of the shear velocity is not unique. Three different definitions will be used. 
The first two definitions are based on the wall shear-stress at each of the two walls, 
that is, 


u T , = u 


dUj 

dr' 


i/2 


U To = 


dUe 

dr 


1/2 


(5.0.1) 


These will collectively be called local u r . The third definition is global; it is obtained 
by analogy to the plane channel. In the plane channel, the mean pressure gradient 


.A •* 
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^ is — 1 when normalized by pu 2 ; in the curved channel, the mean pressure gradient 
is - jy- so the global « T is defined as 


( 1 dP\ l/7 _ (r % 2 u T 2 -|- r o 2 u T 2 0 \ 

V r c pd9j { 4 r? J 


1/2 


(5.0.2) 


The global value u Tg will be referred to simply as u r . The Reynolds number based 
on u T and 8 is 168 for the results presented here. This corresponds to a Reynolds 
number of 2990 based on the centerline mean velocity. Reynolds numbers based on 
u Tl and u TO (Re, and RtJ are 155 and 180 respectively. 

In these computations, 128 Fourier modes are used to represent the velocity field 
in the z and 9 directions. Chebyshev polynomial' through order 64 are used in the 
r direction. In wall units, grid spacing in the z direction is A z + = Azu r ju = 6, 
and in the 6 direction it is r+ A 9 ~ 18. In the r direction, the closest grid point to 
the wall is at = 0.2, and toward the center of the channel the maximum spacing 
is A y + — 8.2. These grid dimensions are for the coarsest collocation grid that 
can be used to represent the Fourier/Chcbyshev velocity representation without 
loss of information; thus the grid dimensions are indicative of the resolution of 
the velocity field (this is not the collocation grid used to compute the nonlinear 
terms). In computing th convective terms, aliasing errors were removed since for 
time dependent problems aliasing may be particularly damaging (Moser, Moin, & 
Leonard 1983). A time-step of 0.00056/u r was used in these computations, which 
yielded a maximum Courant number of 0.8, where the Courant number C is defined 
as 


C = icAt 


( 


vq 

rA9 


Vy 

A r(r) 


+ 


Az 


)■ 


The initial condition for these computations was obtained from a low-Reynolds- 
number, large-eddy simulation of Moin & Kim (1982). The velocity field from their 
calculation was simply interpolated to the collocation grid for the current calcula- 
tion. It was then allowed to evolve for about 12 S/u T , at which time the flow reached 



statistical equilibrium. The calculations were then continued in order to obtain an 
adequate statistical sample. Statistics reported here wore averaged over a time of 
about 6 6/u t , which corresponds to 107 6/U c i, where U c i is the centerline mean veloc- 
ity. The statistics were obtained by averaging in the streamwise direction, in time, 
and often in the spanwise direction. However, averages in the streamwise direction 
and time are not independent because of the Taylor hypothesis. The space time 
correlation functions in §5.5 indicate that the largest eddies arc coherent for a time 
of about 0.5<5 /u t . Thus, the temporal averaging over 66 /u r provides approximately 
12 times better sample of the largest eddies tha i a single velocity field. Also, the 
two point spatial correlations (§5.4) indicate that the largest eddies are coherent 
in the streamwise and spanwise directions over about | the computational domain 
in those directions. The statistics reported here therefore represent approximately 
300 independent samples of the largest eddies. This statistical sample is considered 
marginal for some quantities ( e.g . two point correlations, spectra, and high order 
moments). 

In the sections to follow, we will be concerned with several types of averaging and 
several different velocities. To facilitate discussion, the following notation will be 
used. The velocity vector is denoted v as in the previous sections, with components 
v r , vo, and v z , and ( a) 2 ot denotes a averaged over the z, 6 and t coordinate 
directions. Two special averages are defined as a — ( a) z ot and a = ( a.)ot ■ Several 
averages of the velocity are defined as U, = vi, u, = v, — Z7,, and u' = t\ — 
To facilitate compai son with the plane channel, u, v, and w will be used 
interchangeably with uq, u r , and it z ; for example U - JJo and v' = u' r . Finally, 
the superscript + indicates normalization by local u T and u. All velocities are 
normalized by global u T unless otherwise indicated. When quantities are plotted in 
local wall coordinates, they will be normalized such that positive normal velocity is 
directed away from the wall. 
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5.1 Mean Statistics 


The mean velocity has been plotted in law of the wall coordinates •• ^igure 5.2. 
Both u { and y ' arc based on the local tt T . Also plotted are the piane chanrel d\ta 
of Eckehnann (1971) at Reynolds number Re r = 146. The mean velocities on both 
sides of the channel and the data of Eckehnann are in excellent agreement for y + less 
than 20. For y + greater than 20, the mean velocity of the unstable (concave) side 
lies below the other two. This is the experimentally observed effect of curvature. In 
the experiments of Hunt & Joubert (1979), with approximately the same curvature 
and a 10 times larger Reynolds number than in these .'.omputations, the mean- 
velocity proxies of the concave and convex sides did not diverge until y ’ ft: 200 (see 
Figure 1.1). In both the experiments and the computations, however, the point of 
divergence is at approxiinatly the same y/r location of 0.0015 (herey is distance 
from the wall). Th : s is in accordance with the conjecture of Hoffman & Bradshaw 
(1978) that curvature will not affect the law of the wall until y/r is sufficiently large. 
Also note that the mean- velocity profile of the convex side is in good agreement with 
the plane-channel profile of Eckehnann, indicating that the effect of curvature on 
the mean- velocity profile is not significant on the convex side. In the experiments 
of Hunt & Joubert (1979) the convex wall mean velocity profile was not in good 
agreement with the plane channel profile at large y + . This may be due to the much 
higher Reynolds number of the experiments. 

In strongly carved channels, Wattendorf (1935) and Eskinazi & Yeh (1956) ob- 
served a region of constant mean angular momentum ( rU ). In a region of constant 
mean angular momentum, the mean vorticity is zero (potential fLw in the mean). 
However, in a mild curvature case (6/R « 0.01), Hunt & Joubert (1979) did not 
observe a constant mean angular momentum region. The mean angular momen- 
tum profile from the current calculation (Figure 5.3) also shows no constant region, 
which confirms that it is probably a stiong curvature effect. 
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Turbulent intensities as a function of radial position are shown in Figure 5.4a. As 

expected, turbulence intensities near the concave wall arc substantially higher than 

those near the convex wall. As can be seen in Figure 5.46, this difference persists 

even when the contribution of the Taylor-Gortler vortices if is not included (see 

§5.2). In Figure 5.5 the intensities (u 2 ) are plotted in local wall coordinates for 

both sides of the channel. Also plotted are the plane- channel data of Kreplin & 

Eckclmann (1979a) at Reynolds number Re T = 195. The streamwise intensities 

(u 2 ) for both curved walls and the plane channel are in good agreement when 

normalized in this way. The spanwise (to 2 ) and normal (v 2 ) intensities on both 

curved walls are also in good agreement when normalized by local u r . However, their 

1/2 

agreement with the data of Kreplin and Eckclmann is not as good as that of u 2 
In particular, the computed v intensities are considerably below the experimental 
plane channel data. The reasons foi this are not known, but it is unlikely to be a 
curvature effect. 

The turbulent shear stress (— tiv) is presented in Figure 5.6a along with the con- 
tribution of the Taylor-Gortler vortices to the turbulent stress (-{/{;) and the total 
shear stress (viscous and turbulent). The differences between the concave and con- 
vex sides are striking. In particular, away from the wall the Taylor-Gortler vortices 
make a significant contribution to the concave side Reynolds stress (as much as 
40%), but they contribute negligibly to the convex side. In Figure 5.66 where — u'v' 
normalized by local u r is plotted, it is clear that the curvature has enhanced the 
Reynolds stress on the concave side relative to the convex side. Figure 5.6c, in 
which — uv in local wall coordinates and the data of Eckelmann (1974) are plotted, 
shows that Eckelmann ’s plane channel data at Re T = 146 lies between the concave 
and convex wall Reynolds stress. The correlation coefficient m)/(u 2 v 2 ) 1 / 2 shown 
in Figure 5.7 indicates that streamwise and normal fluctuations arv more corre- 
lated on the concave side than on the convex side (coefficients of 0.5 as opposed to 
0.4). Away from the wall, this is in part a result of the Taylor-Gortler vortices, as 
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is seen in the correlation coeflicient with the contribution of the vortices removed 
(uV/ju^v^) 1 / 2 ). 

On both sides of the channel the correlation coefficient shows a local maximum 
near the wall (y* « 10). The same behavior of the correlation coefficient is dis- 
cernible in the calculations of Moin Sc Kim (1982). This is most probably a result 
of the presence of organized motions near the walls. 

Profiles of rms voiticity fluctuations excluding the contribution of the Taylor- 
Gortier vortices (w' 2 ) normalized by u T /6 are shown in Figure 5.8 a. Because of 

the large spatial scale of the Taylor-Gortler vortices their contribution to the rms 
vorticity fluctuations is negligible (less than 4%). As was observed by Moin Sc Kim 
(1982), the spanwise vorticity profile attains its maximum at the wall, and falls off 
monotonically away from the wall, and the streamwise vorticity profile attains its 
maximum at the wall and has a local maximum at y + zz 20. This streamw'sc rms 
vorticity profile is consistent with the presence of streamwise vortices near the walls 
(see for example Bakewell Sc Lumley 1967; Blackwelder Sc Eckelmann 1979), since 
such streamwise vortices would have vorticity concentrated in their cores, account- 
ing for the maximum in vorticity fluctuations at y + ~ 20. The vorticity would also 
change sign as the wall is approached because of the presence of the wall and the 
no-slip boundary conditions. The change of sign would account for the minimum in 
rms vorticity at y + « 5. The streamwise vortices will be discussed further in §5.4 
and §5.6. A maximum in the streamwise and spanwise rms vorticity at the wall is 
also an expected consequence of the “splatting” effect (Moin Sc Kim 1982) to be 
discussed in §5.3. Away from the walls, the three components of the rms vorticity 
fluctuations are virtually identical, in contrast to the velocity fluctuations, which 
axe significantly different. As explained by Moin Sc Kim (1982), the contribution of 
small-scale fluctuations to the rms vorticity is much larger than their contribution 
to the intensities, and away from the wal. the small-scale fluctuations are expected 
to be isotropic. 


45 


In Figure 5.86 the rms streamwise vorticity fluctuations, nondimensionalized by 
local u*/u for each wall, are shown, as well as the plane-channel data of Kastrinakis 
& Eckelmann (1983), at Re r — 580. The profiles from both walls are in very good 
agreement when nondimensionalized in this way, and they arc in good agreement 
with the plane-channel data for y * greater than 10. However, the experimental 
plane-channel profile does not obtain a minimum near the wall, and the computa- 
tional and experimental limiting wall values are in disagreement. Other researchers 
have observed a range of limiting wall values of w x from 0.065 to 0.12 (see Kreplin 
& Eckelmann 1979a) and Moin & Kim (1982) calculated a value of 0.13; 0.19 was 
calculated here. The reason for this discrepancy is not known, but it is unlikely to 
be a curvature effect since it is the same on both curved walls. The plain-channel 
calculations of Moin & Kim (1982) and the present calculations, which use un- 
related numerical methods, both show local minima in rms streamwise vorticity 
fluctuations near the wall. This suggests to us that the computed results, which 
show the minimum, may be correct despite their disagreement with Kastrinakis 
and Eckelmann. However, it is possible that the results of the current calculations 
axe affected by the projection error (see §2.2). This uncertainty precludes a final 
conclusion on the near-wall behavior of the streamwise vorticity. The limiting wall 
value of the spanwise vorticity fluctuations (0.36) in the present calculations is also 
higher than observed experimentally (0.2 to 0.3, Kreplin & Eckelmann 1979a) and 
the computed value (0.2) of Moin & Kim (1982). 

In Figure 5.9 the skewness factors of it, v, and it; (S( it) = u 3 /u 2 ) are shown 

with and without the contribution of the Taylor-Gortlcr vortices; (S(u) and S(u')). 
When the Taylor-Gortler vortices are included, the skewness of it becomes large 
and negative away from the walls (as large as -77.22). The reason for this behavior 
will be discussed in the following section. Because of the reflection symmetry of the 
Navier-Stokes equations, the skewness of w should be zero. The very small values 
of to skewness shown in Figure 5.9 indicate that the statistical sample from which 
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the skewness is calculated is adequate. In Figure 5.10, skewness factors from both 
sides of the channel are plotted in local y* coordinates, together with the data of 
Krcplin & Eckelmann (1979a) for the plane channel. Note that the u' skewness 
factors are in very good agreement with the plane-channel data. The agreement for 
the v' skewness is not nearly as good. The v' skewness of Kreplin and Eckelmann 
never becomes negative and has a much larger value at the wall. Also, recent S(v) 
measurements by Alfredson & Johansson (1984), which were limited to y + > 30, 
show no tendency to become negative near the wall. Moin & Kim (1982) also 
observed negative v skewness factors in the vicinity of the wall. 

It is interesting that the it' skewness at the wall is approximatly 60% higher on 
the convex side than on the concave side, indicating that the large u fluctuations 
associated with high speed fluid arriving from away from the wall are stronger on 
the convex side. This may be attributed to the effect of the Taylor-Gortler vortices 
on the underlying turbulence. On the concave, side there is a region of strong flow 
away form the wall (see §5.2), which would tend to inhibit the motion of high speed 
fluid toward the wall. 

Velocity flatness factors of u, v, and w (F(u) = u 4 /u 2 ) are shown, with and 
without the contributions of the Taylor-Gortler vortices (F(u) and F(u’)), in Figure 
5.11. With Taylor-Gortler vortices included, the flatness of u is very large away 
from the wall (see §5.2). When the contribution of the vortices is removed, the 
flatnesses of all three velocity components are between three and four away from 
the walls (a Gaussian distribution has a flatness of 3). Near the wall, the flatness 
factors generally become large, which is indicative of int *rmittence or spottiness 
of turbulence near the walls. The u' flatness factors for both curved walls are in 
good agreement with the plane-channel data of Kreplin and Eckelmann (Figure 
5.12). The w' and particularly the v' flatness factors do not agree as well with the 
experiments. The computations of Moin & Kim (1982) show similar disagreement 
of v' flatness with the data of Kreplin and Eckelmann. Near the walls, the flatness 
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of v' and to' arje extremely large (as high as 30.83 for v 1 and 9.62 for w'). Also, the 
flatness factors of all three velocity components are higher at the convex wall than 
at the concave wall, which suggests that the turbulence is more intermittent very 
near the convex wall. This may be a consequence of the lower Re r on the convex 
side. 

In order to evaluate Bradshaw’s (1973) suggestion for a correction to the mixing 
length for curved flows, the constant /? in the expression, 


l = 1 .p 

<0 H dU/dr - U/r 

was calculated. Here, t is the mixing length computed from 

1/2 


(5.1.1) 


1 = 
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l i (W _ U\ 

r V dr r / 


(5.1.2) 


This definition of mixing length is not meaningful where either tit; or dU/dr — U/r 
is zero (y = -0.163 and 0.033 respectively) because they do not vanish at the same 
place as required by the mixing-length assumption. The reference length to was 
assumed to be that of a plane-channel computed from the expression 


_ * [ _ / _ (1 - \y\) \ 

6 3 V 6 ) 


-e* + / A+ ), 


(5.1.3) 


which is an adaptation of a length-scale expression used by Norris & Reynolds 
(1975). The constants k and A + were chosen to fit the mean-velocity profile of 
Eckelmann (1974). Thus, this reference length can be viewed as the mixing length 
implied by Eckelmann’s mean-velocity profile. However, near the wall where viscous 
stresses dominate the turbulent shear-stress, the mean-velocity profile is not sensi- 
tive to the mixing length, and the reference length lo is not expected to be relevant 
there. The computed values of (3 for the convex and concave walls are shown in 
Figure 5.13. In the region where we expect (5.1.l)-(5.1.3) to be meaningful (say 
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40 < y + < 120), (i varies from about 20 to 2.5, and is consistently lower on the 
convex side than on the concave side. This is consistent with the recommendations 
of Bradshaw (1973). 

It is interesting that on the convex side the value of /? is large despite the fact that 
the velocity profile on the convex side is in good agreement with that of Eckelmann 
(1974) see Figure 5.2. This is because the turbulent shear stresses (— uv) for the 
convex wall and for the plane channel differ by as much as 20%. These shear stresses 
result from virtually identical mean- velocity profiles because of differences in the 
corresponding mean-velocity equations, which result in different expressions for the 
turbulent shear stress in terms of the mean velocity. For the plane channel, the 
shear stress is 


_ . ,dP i du 

where y w is distance from the wall. For the curved channel, 


(5.1.4) 


_ ISP 

UV = 

2 80 
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2 r d U/ r 

Re dr 


(5.1.5) 


where velocities and the pressure have been normalized by the u T local to the convex 
wall. Thus, curvature has a significant effect on the mean equation. 


5.2 Taylor-Gortler Vortices 

In laboratory experiments, Taylor-Gortler vortices can be made stationary by 
introducing weak disturbances into the boundary layer upstream of the curved 
section (see §1.1). These controlled disturbances have the effect of triggering th ’ 
Gortler instability, causing the vortices to grow in preferred locations. In the present 
computation the analogous upstream disturbances are the Taylor-Gortler vortices 
themselves as they are convected out the downstream end of the computational 
domain and are reintroduced at the upstream boundary by the periodic boundary 
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conditions. A similar phenomenon occurs in high-tteynolds-number Taylor-Couette 
flow; although the flow is fully turbulent, there are stationary, axisymnietric Taylor 
vortices present (Coles 1965) because of the periodicity in the azimuthal direction. 
Note that nothing precludes the vortices from moving in the spanwise direction; 
the spanwise periodicity does not force them to remain stationary. However, the 
periodic boundary conditions in the spanwise direction do have the artificial effect 
of restricting the possible wavelengths of the Taylor- G or tier vortices. The results in 
this section concerning the effects of presumably stationary Taylor-Gortler vortices 
are expected to be valid for nonstationary vortices, as long as they are coherent over 
distances and times much larger than the length and time scales of the underlying 
turbulence. Also note that these calculations are fundamentally different from the 
experiments in which disturbances are introduced to lock in the vortices, because 
here no artificial disturbances were introduced. The computations were started with 
a turbulent velocity field taken from the computations of Moin & Kim (1982) which 
was allowed to evolve in the curved channel. The Taylor-Gortler vortices in this 
computation developed from turbulent fluctuations with a broad spectrum. 

In order to study turbulent Taylor-Gortler vortices they must be differentiated 
from the underlying turbulence. For this study the vortices are determined to be 
the average of the velocity field in 9 and t minus the average in. 9, z and 1 (v — v). 
Other definitions axe possible; for example, any spatial and/or temporal filter could 
be used (the averages used here are a special case). Note that since the temporal 
average is over a finite time, the vortices that survive this averaging may actually 
be moving or evolving on a time scale of the averaging time (6 S/u T ) or longer. 
Therefore, with the current method it is not possible to determine whether the 
vortices are drifting or not. If slowly drifting Taylor-Gortler vortices were present, 
the results of the averaging performed here would underestimate their strength and 
effects. Thus the effects of the vortices reported in this and other sections may be 
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understated. Qther schemes for differentiating the Taylor- Cortlcr vortices from the 
underlying turbulence are a topic for future work. 

The Taylor-Cortler vortices were isolated by determining the average velocity v as 
a lunction of r and z. In Figure 5.14 the secondary flow streamlines of the vortices 
are plotted. The streamlines before and after averriging over the mirror image 
flow arc shown ( z is mapped to —z, see §3.3). Note that the effect of averaging 
over the mirror image is to make the contour lines somewhat smoother, and to 
remove a minor asymmetry of the vortices. Averaging over the mirror image flow 
was performed for all the remaining results in this section, which had the effect of 
removing similar asymmetries from the results. In this and all subsequent contour 
plots, negative quantities are denoted by dashed lines. The streamlines show that 
the vortices fill the entire channel, though they are concentrated somewhat on the 
concave side. Between the two vortices, where the streamlines are closely packed, 
is a region of relatively strong flow away from the concave wall. Note that the flow 
toward the concave wall due to the Taylor- Gor tier vortices is significantly more 
diffuse than the flow away from the concave wall. The contorted shape of the 
outer contours is an artifact of the finite statistical sample used; so no particular 
significance should be attached to it. 

In Figure 5.15 the spanwise variation of the wall shear-stress is shown for both 
walls (j^ a (^+ u . ) i ). The z position of the plot is aligned with Figure 5.146. 
On the concave wall, there is a very sharp minimum in shear stress located between 
the vortices. The oscillatory behavior of the shear-stress curves on both walls is 
attributed to a poor statistical sample. On the convex side, the effect of the vortices 
is so small that it is masked by the statistical noise. 

Contours of the average velocities associated with the Taylor-Gortler vortices are 
shown in Figure 5.16. The intense region of negative r velocity is evident in the area 
between the vortices. Note that the largest radial velocity is 0.85u r or about 5% of 
the centerline velocity. This strong radial flow convects low-speed fluid away from 



the wall, giving rise to an area of strongly negative uo that has a magnitude as large 
as 2.8 u T (about 15% of the centerline velocity). It is clear that the region of strong 
negative v is responsible for most of the Reynolds stress contributed by the Taylor- 
Gortler vortices, as can be seen in Figure 5.17 where the contours of tit; are plotted. 
In the middle of the region of strong radial flow, the local Taylor-Gortler Reynolds 
stress is as high as 1.5ttJ (recall that the maximum contribution of the vortices to 
Reynolds stress is about 0.2u 2 ). Also of interest is the significant gradient of the 
streamwise and spanwisc velocities (du/dr and dw/dr ) near the concave wall, as 
indicated by the concentration of contour lines in Figures 5.16a and c. The gradient 
of the streamwise velocity is responsible for the large defect in shear stress on the 
concave side (Figure 5.15). 

It is now clear why the skewness S(u) and flatness F(u) that include the contri- 
bution of the Taylor-Gortler vortices are so large away from the walls. The intense 
region of negative u makes the velocity distribution of the vortices extremely skew. 
This region also contributes greatly to ti 4 , resulting in a very large flatness. There 
is probably no significant effect on S(t>) and F(v) because v is small compared to 
the fluctuations of the underlying turbulence. 

The Taylor-Gortler vortices affect the underlying turbulence by convecting it 
along the streamlines in Figure 5.14, and by introducing a secondary strain field 
which contributes to its production. In Figure 5.18 the contribution of the un- 
derlying turbulence to the components of the Reynolds stress tensor are shown as 
a function of r and z. Plotted are contours of — uju' ; the mean value is 
subtracted to make the variations more apparent. 

In the plots of the diagonal elements of the Reynolds-strcss tensor (Figure 5.18 

- ■ i — «*-**«/ 

a, b and c for -u'g — u'q, u — tzj. 2 , and u'j 2 — u'£, respectively) there is a strong 
positive region slightly away from the concave wall centered on the region of strong 
negative u r (labeled A in the figures). This is the result of the turbulence near the 
wall, where the intensities are maximum (see Figure 5.46), being convected toward 
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the center of the channel by the strong u r . The positive region of u'q — u q is 
the strongest of the three, since the maximum in intensity near the wall is most 
pronounced for the 0 intensity. Likewise, u' 2 — u' r 2 has the weakest such positive 
region, because the maximum in r intensity is least pronounced. Toward the sides 
of the plot domain, where u r is weak but positive, the opposite mechanism (fluid 
with a low turbulence level convected toward the concave wall) produces the regions 
of negative u' 2 — it' 2 in the region labeled B. It could also be argued that a similar 
convection mechanism is at work near the convex wall; however, the effect is much 
weaker and cannot be reliably differentiated from statistical noise. 

Near the concave wall (y 4 < 20) there is a region of very intense negative u'q -u'q 
under the positive regions discussed in the previous paragraph (labeled C). There 
is a similar, though considerably weaker, region of negative u'j 2 - u'}. In Figure 
5.16a we saw that duofdr was positive in this same area. Thus, production of u'q is 
adversely affected, which would contribute to the depression of u'£ -tx' z 2 in region C. 

The Reynolds shear-stress term u' 0 u' r — u' 0 u' r (Figure 5.19) is similar to u'q - u'q 
in that it is positive in region A and negative in regions B and C. Convection is 
responsible for the behavior in regions A and B. In region C, the production of u' 0 v' r 
is suppressed, contributing to the negative values there. 

5.3 Reynolds Stress Balance Equations 

The Reynolds-stress equations in cylindrical coordinates are derived in Appendix 
C. Here we consider the Reynolds-stress equations for the special case in which 
the mean velocity Uq only varies in the radial direction. For this special case, the 
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The equation for the turbulent, kinetic energy q 2 = ^(u 2 + v 2 + u> 2 ) is 


dq 2 __ dU/r 

3t ~ UVT dr 


1 drvq 2 — =- 

°-2r-ar ~ U VP 


+ JL 1 1 

21te \ r dr dr 


(5.3.2) 


where e is the dissipation of turbulent kinetic energy. In these equations ((5.3.1) 
and (5.3.2)) the terms on the right-hand sides are labeled (in order of appearance) 
production, convection, turbulent diffusion, velocity-pressure gradient (VPG ), viscous 
diffusion, and dissipation. Zeros appearing in the equations indicate terms that are 
identically zero. Many of the terms in the Reynolds-stress equations in cylindrical 
coordinates do not appear in the corresponding equations in Cartesian coordinates; 
they are the so-called “extra” terms (Bradshaw 1973). These terms reflect the fact 
that in cylindrical coordinates the orientation of the coordinate axes is a function of 
6. If the flow is homogeneous in the 0 direction, the orientation of the mean velocity 
vector and the principal axes of the Reynolds-stress tensor are also a function of 9. 
This gives rise to streamwise (0) gradients of mean velocity and Reynolds stress. 

The streamwise gradients contribute to production, convection, and diffusion of 
the Reynolds stresses. For example, in the it 2 equation the production term consists 
of two parts: -2uvr~£- which represents the production of u 2 by interaction of 
turbulence with the mean shear; and ~2uvj which is the production caused by the 
interaction of turbulence with the streamwise gradient of the mean velocity vector. 
Similar streamwise production terms appear in the balance of v 2 and uv. The 
convection terms in each of the equations represent the convection of the Reynolds 
stresses by the mean ’.eocity; this is not zero because of the nonzero streamwise 
gradients of the Reynolds-stress tensor. The diffusion terms consist of diffusion in 
the radial direction and diffusion in the streamwise direction, which acts to 

diminish the gradient of the stress tensor in the streamwise direction. 
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In the ti 2 equation, all the streamwise gradient terms are the same as in the v 2 
equation with opposite sign, and there arc no streamwise gradient t' 1 . in the w 2 
equation. Thus, these terms do not contribute to the equation for the kinetic energy 
q 2 . The reason is that there is no stvcainwise gradient of q 2 since it is an invariant 
of the Rcynolds-stress tensor. 

Equations 5.3.1 are derived from the Navier-Stokes equations. When the Navier- 
Stokes equations are discretized in space and time by the method outlined in §2 and 
§3. an additional term enters each equation, which represents the contribution of 
the projection error (see §2.2). This projection-error term, U{P(A)j -f ujP(A)i (A 
is as defined in (2.1.12)), was computed and found to be uegligibly small compared 
with other terms in the equations. 

In Figures 5.20, 5.21, 5.22 and 5.23 the various terms in the Reynolds-stress bal- 
ance equations (5.3.1) are plotted in local wall coordinates for both walls (velocities 
nondimensionalized by local u r and lengths by i//u r ). This nondimensionalizatio. 
is consistent with the wall-similarity hypothesis, and attempts to eliminate the ef- 
fect of the different Reynolds numbers on the concave and convex walls. Except 
for the streamwise turbulent diffusion in the uv equation, none of the terms due to 
streamwise gradients are included because they are negligibly small. Terms in each 
of the equations show remarkably little difference between the concave and con- 
vex sides of the channel when plotted in local wall coordinates. The few significant 
differences will be discussed after we examine the common features. 

The u 2 equation is largely dominated by production and dissipation. There is 
a large peak in production near the wall ( y * « 15), which is balanced in part by 
the large dissipation near the wall. Turbulent and viscous diffusion carry u 2 energy 
from the region of maximum production (note the minima in viscous and t urbulent 
diffusion) in both directions, away from and toward the wall. Very near the wall, 
the extreme values of dissipation are balanced by diffusion from the maximum 
production region. Far from the wall, production and the positive contribution of 
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turbulent diffusion are balanced by dissipation and the velocity pressure gradient 
terms. I r the u 2 equation, the velocity pressure gradient term consists entirely 
of the pressure-strain correlation which represents transfer of energy to the other 
components of turbulent intensity (v 2 and in 2 ). 

Since the v 2 and w 2 equations contain no significant production terms, their only 
source of energy is the pressure-strain correlation. In Figure 5.24, the pressure-strain 
terms appearing in the u 2 , v 2 , and w 2 equations are plotted together. Beyond y + 
of 20, the major effect of the pressure-strain correlation is to distribute energy from 
the tz 2 component to the v 2 and tv 2 components. However, close to the wall there is 
a large transfer from the normal component, v 2 , to the other compos its. This was 
observed by Moin & Kim (1982) in their computations of plane-chr nel flow and was 
referred to as the “splatting” or impingement effect. It is caused by fluid elements 
coming from away from the wall, impinging on it, and transferring their energy to 
motion parallel to the wall. Because of the no-slip boundary conditions and vortex 
stretching, the splatting effect gives rise to large streamwise and spanwise vorticity 
fluctuations, as seen in Figure 5.8. 

In the v 2 equation, the pressure- strain and pressure-diffusion terms combir" to 
form the velocity pressure gradient term. In Figure 5.21 the velocity pressure- 
gradient correlation is the major positive contribution to v 2 . Note that it is only 
slightly negative near the wall, implying that the pressure diffusion term is posi- 
tive near the wall to cancel the negative pressure-strain term. The pressure-strain 
correlation, which is the source of v 2 energy, is maximum at y + ss 35. As in the 
u 2 equation, energy is diffused from this location in both directions, toward and 
away from the wall, the predominant term being the turbulent diffusion. Kinematic 
constraints on the normal velocity ( dv/dr — 0 at the walls) require that the vis- 
cous diffusion and viscous dissipation of v 2 have zero slope at the wall. This is not 
apparent in Figure 5.21; however, when the region around the origin is magnified 
(Figure 5.25) it can be seen that these slopes are indeed zero. 
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Pressure-strain and viscous dissipation dominate the w 2 equation. Very near the 
wall, however, there is significant viscous diffusion. 

In the uv equation, the production dominates, and again there is viscous and 
turbulent diffusion from the maximum source region (y f ~ 15) toward and away 
from the wall. However, in this case the viscous dissipation is negligible almost 
everywhere, and the production is balanced by the velocity pressure gradient and 
turbulent-diffusion terms. 

It is interesting to compare the current Reynolds-stress balances to those obtained 
by Moin & Kim (1982). As with other statistical correlations, there is a remarkable 
similarity of the u 2 and w 2 balances in this study and those reported by Moin and 
Kim, though there are differences in the y + locations of the maxima, minima, and 
zero crossings of the various terms. In general, the y f location of each feature is 
larger in the calculations of Morn and Kim. This difference is due to the resolution 
in their computations which was inadequate to resolve the wall-layer structures at 
their proper scale. 

The v 2 and uv balances appear quite different in the two calculations. However, 
in the case of v 2 , if we approach the wall from the center of the channel, the 
same features a observed in both calculations though at different y x locations. 
In addition, in the vicinity of the wall ( y + < 15), Moin and Kim show a relatively 
large magnitude of turbulent diffusion balanced by a large velocity pressure-gradient 
term. This is not found in the present calculations. Note that in both calculations 
the location of the maximum in the turbulent diffusion term ( y + = 15 here and 
y v = 30 in Moin and Kim) is at the same location as the minimum in the v 
skewness factor. 

The balance of the turbulent kinetic energy ‘zq 2 = u 2 -\-v 2 + w 2 is shown in Figure 
5.26. The kinetic energy equation is dominated by the tt 2 term, so this balance is 
very similar to the u 2 balance. As w?" seen above the turbuk. t diffusion is positive 
very near the wall as a result of t,:v i'usion of energy from the maximum source 



region. In contrast,, the estimated balance of Townsend (1970) shows no positive 
region of turbulent diffusion. Moreover, Townsend shows a very large pressure 
diffusion term near the wall which is also contrary (in relative magnitude) to the 
current results and to those of Moin &; Kim (1982). Townsend’s estimates for the 
remaining terms arc in qualitative agreement with the current calculations. 

As was noted above the terms of the Reynolds-stress balance are remarkably simi- 
lar on the convex and concave sides (when normalized by local wall variables). There 
are, however, several significant differences. In the u 2 equation, the production is 
somewhat higher on the concave side (about 10%), and near the wall (y + < 25) the 
turbulent diffusion is lower. The opposite is true in the v 2 balance, in which the 
turbulent diffusion is larger on the concave side. The uv balances show the most 
differences between the concave and convex sides. This is not surprising since uv 
itself shows more differences between the two sides than the turbulence intensities. 
On the concave side the peak production of -uv is about 5% less than on the convex 
side, a result of the smaller values of dU/dr on the concave side when expressed 
in local wall coordinates. The velocity pressure gradient term is as much as 20% 
greater on the convex side, and the turbulent diffusion from the maximum source 
region is about 40% lower on the concave side. Also, in this balance the streamwL i 
turbulence diffusion is significant and contributes to — uv on the concave side and 
diminishes -uv on the convex side. Note that away from the walls the streamwise 
turbulent diffusion is as large as the radial diffusion (about 20% of production). 

Many of the. differences cited above are in the turbulent diffusion and pressure 
strain terms. The turbulent diffusion terms in these calculations include several 
effects: the convection of the underlying turbulence by the Taylor-Gortlcr vortices, 
the actual turbulent diffusion of the underlying turbulence, and the enhancement (or 
diminishment) of that diffusion by the Rryleigh instability mechanism. The effects 
of the Taylor-Gortler vortices and the Rayleigh mechanism on the concave wall 
will be opposite that on the convex wall, so it is not surprising that the turbulent 


diffusion is different on the two walls. This is in accordance with Bradshaw’s (1973) 
assertion that curvature effects on the lleynolds-stress equations must appear in the 
higher order statistical correlations. Bradshaw’s argument would also suggest that 
the pressure strain terms should be affected since they can be expressed as integrals 
of two-point triple correlations of the velocity gradients (see for example Launder, 
Reece, & Rodi 1974). The significant curvature efTects on the uv pressure-strain 
terms is also in accordance with the suggestions of L.aundcr, Reece, &; Rodi (1974) 
and So (1975) that curvature effects can be accounted for by properly modeling the 
pressure strain terms. 

It is interesting that the dissipation terms, which are dominated by the small 
scales, are in very good agreement on the two walls. Also, as noted in §5.1, the 
rms vorticity fluctuations, which are sensitive to the small scales, were in good 
agreement in local wall variables. This suggests that curvature has a minimal effect 
on the small scales of turbulence. 


5.4 Spectra and Two Point Spatial Correlations 


One dimensional, two-point correlation functions, 




- =3 


(5.4.1) 


(no summation) are plotted in Figures 5.27 and 5.28 at six r locations. The corre- 
lations have nonzero values for large separations, indicating that the length of the 
computational domain in the streamwise direction is somewhat inadequate. How- 
ever, the values are small, which suggests that the periodic boundary conditions in 
the streamwise and spanwise directions have not significantly distorted the compu- 
tations. Also, as expected, for small displacements, the correlation of the velocity 
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component in the direction of the displacement is greater than for the transverse 
components. 

Very near the walls, Roo(68) decays very slowly with increasing 66, which indicates 
that structures are coherent over long distances in the strcamwisc direction. In the 
spanwise direction, Rgo(Sz) decays to zero very quickly, showing that the spanwise 
extent of these structures is much smaller. Also, Roo{6z ) becomes negative and 
reaches a minimum at 6z + as 50. The presence of a negative minimum in R#o(6z) 
and the slow decay of Rgg(66) are clear evidence of elongated regions of high and 
low speed fluid alternating in the spanwise direction (streaks). The location of the 
minimum is the mean distance between a high speed and a low speed streak; thus, 
the mean streak spacing is twice this distance (A + as 100), in agreement with Kline 
et al. (1967). At a greater distance from the walls, the pronounced minimum in 
Roo{6z) does not appear. This is also in accordance with the observations of Kline 
et al. that the high and low speed streaks are confined to a region close to the wall. 

Note that near the wall, R rr (6z) and R zz (6z) also attain negative minima, indicat- 
ing the existence of alternating regions of positive and negative u' r and u' z , though 
the streamwise coherence of these regions is significantly smaller than that for u' e . 
The location of the minimum in R zz (6z) is at 6z rs 50, but the minimum in R rr [6z ) 
is at 6z « 25. 

Bakewell Sc Lumley (1967) and Blackwelder Sc Eckelmann (. ‘'"'9), among others, 
have proposed that the wall region is dominated by pairs of com. rotating stream- 
wise vortices. These vortex pairs are thought to have a long streamwise extent, and 
to be responsible for the low and high speed streaks. However, the streamwise two- 
point correlations show that the streamwise coherence of v' and w' fluctuations is 
smaller than that of u' fluctuations; as will be seen in §5.6, ths v' and w' velocities 
do not exhibit streakiness near the walls. If such vortex pairs provided a significant 
contribution to v' and w' fluctuations, they would give rise to spanwise two-point 
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correlation functions of v and w with negative minima. The position of the mini- 
mum in R zz (8z) would be related to the mean distance between vortices of opposite 
sense, and the minimum in R rr (6z ) would be related to the mean distance across 
a vortex. In Figure 5.29 the spanwisc two-point correlation functions at several y 
locations near the concave wall are shown. An important feature of these correla- 
tion functions is that R zz {6z) does not have a negative minimum beyond y ¥ — 18, 
whereas Roq(6z) and R rr {8z) do. The absence of a minimum in R zz (6z) implies that 
vortex pairs, with their centers at the same y location, are not a dominant feature of 
the near-wall turbulence. But, this does not preclude the existence of vortex pairs 
with the vortex centers at different y locations. The existence of solitary vortices 
is consistent with the behavior of R rr (Sz). Moreover, the fact that the location of 
the minimum in R rr {8z) moves to larger separations as the distance from the wall 
increases, indicates that vortices away from the wall have larger diameters than 
those near the wall. 

The spanwise two-point correlations reported here are in disagreement with the 
measurements of Gupta, Laufer, h Kaplan (1971) in a turbulent boundary layer. 
Their experimentally determined R xx (8z) measured at y + < 12 did not become 
significantly negative and did not have a distinct minimum. They showed no indi- 
cation of the presence of the near- wall streaks. This led them to abandon R xx (6z) 
as a useful quantity for investigating the near-wall streaks. However, it is clear that 
if the streaks are the dominant structures in the viscous sublayer, as observed by 
many investigators, they must produce spanwise two-point correlation functions of 
u' with a negative minimum as presented here. These calculations, the experiments 
of Bakewell & Lumley (1967) and Tritton (1967), and the computations of Moin 
& Kim (1982) all obtain minima in the spanwise two-point correlation of u\ Thus 
we must conclude that the correlation functions of Gupta et al. and their conclu- 
sion that the two-point correlation function carries no structural information tire in 
error. 
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Length scales related to the integral length scale may be determined from the 
one-dimensional, two-point correlation functions using the following definitions: 


/to(» = 2£?/r) = 0.1 , 
Ru(6z = 2 L?) = 0.1 . 


(5.4.2) 


The choice of 0.1 is arbitrary. This definition is used instead of the integral scale 
because it is insensitive to the correlation function at large separation, where statis- 
tical noise is significant. The values of these length scales are plotted in Figure 5.30 
along with a length scale used by Norris & Reynolds (1975) for the plane channel. 
Note that the strearuwise length scales turn up near the walls (especially L ° 0 ), and 
that the spanwise length scales are significantly smaller near the walls than in the 
regions away from the walls. Also, L 0 and £* are greater than the other length- 
scale components in accordance with the observation that correlations are greater 
for components in the direction of the displacement. It is also of interest to calculate 
L\ to compare with the streamwise and spanwise length scales. 

The two-dimensional, two-point correlation functions, 


Ru{» t Sz) 


u[ (r, 9, z, f)u(- (r, 9 + <50, z + 6z, t) 


u 


1 2 


(5.4.3) 


have been computed at selected r locations. Isocorrelation contours are shown 
in Figures 5.31, 5.32 and 5.33. Near the walls, these correlations are elongated 
in the streamwise direction, while surrounded by negative regions in the spanwise 
direction. The positive region of the Rgg correlation extends approximately 400 wall 
units in the streamwise direction, which is indicative of the length over which the 
high and low speed wall streaks are substantially coherent. Note that the R rr and 
R zz correlations arc not as coherent in the streamwise direction. Also, the width 
in the spanwise direction of the R rr correlation is considerably smaller than the 
other two. These observations are consistent with the deductions made from the 
one-dimensional correlations discussed above. 
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Away from the wall, only the Roo correlation is significantly elongated in the 
streamwise direction; R rr and R zz have contour lines which arc almost circular. 
This is consistent with the spanwise space-time correlations of Kovasznay, Kibens, 
ft Blackwelder (1970), which are related to the space-space correlations through the 
Taylor hypothesis. 

One-dimensional energy spectra, defined as the Fourier transform of the two-point 
correlations, have been computed (they appear in Figures 5.34 and 5.35): 



u'(r,0,z,t)u'(r,0 + 9',z,t)e ike °' dO' 
u'(r, 0, z, t)u'(7-,0,z + z',t)e ,k,z ' dz' 


I 


(5.4.4) 


All of these spectra show at least a 2 decade reduction of the high wave-number 
energy density compared with the low wave-numbers. In addition, none of the spec- 
tra show a significant upturning at the highest wave-numbers. These observations 
indicate that the present computational resolution is adequate. 

Throughout the channel, the spectra at high wave numbers (as well as the low 
wave numbers) exhibit anisotropy. Near the walls, the spectra E^k^) attain local 
maxima. These maxima occur at the wave numbers corresponding to the minima in 
the two-point correlation function ( k z « 10 for E$o and E zz , and k z as 20 for E rr ; 
here we are discounting the small glitch in E zz (k z ) very near the concave wall as 
statistical noise). The spectra E rr (k z ) exhibit a local maximum throughout most of 
the channel, just as R rr (6z) maintains its negative minimum far from the wall. Far 
from the wall, the spectra Eu(ko) are not smooth for low wave numbers, which is an 
indication of poor statistical sample. This is probably a result of the slow evolution 
of large structures away from the walls, and may be related to the evolution of tae 
Taylor- Gor tier vortices. 

The dissipation spectra are shown in Figures 5.36 and 5.37. The dissipation spec- 
tra are defined such that their integrals over k are the dissipation given in (5.3.1). 



In these spectra, the values of the maximum dissipation at low wave numbers are 
a decade larger than the values at the high wave numbers. The spanwise spectra, 
however, generally exhibit a substantial upturn at high wave numbers, indicating 
some deficiency in the spanwise resolution. The amount of dissipation represented 
in these upturned tails, which is a measure of the dissipation in the unresolved 
small scales, is generally less than 10% of the total dissip xtion. Thus, most of the 
viscous dissipatiou resides in resolved scales, which is an indication of the overall 
adequacy of the computational resolution. Note, however, that we do not resolve 
the Kolmogorov length scale (rj + « (Re/C/ ci ) 1/4 « 2) in this calculation. Most of 
the viscous dissipation takes place at scales larger than 15rp Also of interest is the 
anisotropy of uhe dissipation at high wave numbers, similar to the anisotropy in 
energy; this may be due to the low Reynolds number of this computation. 

Spectra in the normal direction can be defined from the Chebyshev representation 
of the velocity field, 

N 

My> d >z,t) = J 2 a i{n,0,z,t)T n {y ). (5.4.5) 

n=0 

The Chebyshev coefficient spectra are defined as 

C,(n) - a?(n). (5.4.6) 

These spectra (Figure 5.38) show a 3 to 4 decade fail off from low to high n. However, 
C z has a substantial 1 ’pturn at high n, indicating a difficulty in the resolution of 
u x . It is suspected that this difficulty is related to the projection error (see §2.2). 

5.5 Temporal Spectra and Correlations 

The temporal energy spectra and space-time correlations have been obtained 
using the method described in §3.3. In Figure 5.39 one-dimensional temporal energy 
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spectra are plotted at four r locations. Shown are the spectra determined at the 
chosen frequencies (see §3.3). Note that the spectra start to turn up at the two 
highest frequency points, the reason being that the values of the spectra at these 
points are of the same order as round-off errors. Kourd-off error is one part in 10 7 
because the data were packed to 32- bit words to save storage. These last points have 
not been included in the interpolations used to obtain the two-point correlations. It 
is also interesting that the spectra decay very rapidly beyond J? « 300 (see below). 

The temporal spectra, rescaled using the Taylor hypothesis based on the local 
mean velocity, are 

K (K = ~n) = , (5.5.1) 

where Lt and Lo are the lengths of the domain in which the transforms are taken 
in time and 6, respectively; the spectra are shown in Figure 5.40 along with the 
streamwise spectra at the same r locations. The local mean velocity was used for this 
purpose instead of any of the convection velocities discussed below, because it has 
been observed that °mall scales are converted at about the mean velocity (Sternberg 
1967). Thus, the use of the mean velocity should bring the high wave-number 
part of the temporal and spacial spectra into agreement. Convection velocities 
appropriate for the low wave-number part of the spectra were not used because 
statistical noise precludes accurate comparisons of the low wave-number spectra. 
At all the r locations, the spectra at high wave numbers are in remarkably good 
agreement, indicating that the small-scale turbulence is indeed convected with the 
local mean velocity. At the two r locations nearest the walls (y + — 26 and 37), the 
low wave-number parts of the spectra are also in good agreement. Farther from the 
wall, the low wave numbers appear to be only in fair agreement, though the issue is 
clouded somewhat by statistical noise. This will be discussed further below. Note 
that the location of the downturn in the temporal spectra mentioned above is at the 
cutoff of the computed streamwise spectra when scaled using the Taylor hypothesis. 
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Beyond this wcive number, the computed streamwise spectra are formally zero; thus, 
by the Taylor hypothesis, the temporal spectra should also be zero. This results in 
the sharp downturn in the temporal spectra. 

Contours of the streamwise space-time correlation functions at several r locations 
are shown in Figures 5.41, 5.42, 5.43 and 5.44. Also drawn in these plots is a line 
with slope equal to the local mean velocity. If the Taylor hypothesis were satisfied 
exactly, all isocorrelation contours would be straight parallel lines, and their slope 
would be the convection velocity. In the contour plots we see that near the origin the 
contours are nearly parallel to the line representing the mean velocity. However, 
far from the origin the elongated region of positive correlation bends below the 
mean-velocity line, indicating that structures which are coherent over longer times 
are convected at a reduced velocity. It has been found by several authors (Favre, 
Gaviglio, &; Dumas 1957; Willmarth 1975; and Wills 1964) that convection speed 
varies with the size of the structure being convected. Sternberg (1967), using the 
boundary-layer data of Favre et ai, found that convection speeds for large-scale 
disturbances are as much as 25% greater than the local mean velocity very near 
the wall, and about 25% lower than the mean velocity far from the wall. Small- 
scale structures were found to convect with the mean velocity. Since structures 
that are coherent over long times are expected to be large, it is not surprising that 
our contours indicate that their convection velocities are different from the mean 
velocity. 

A convection velocity may be defined as the slope of the curved line which traces 
the path of slowest descent of the correlation function (Willmarth 1975); in this case 
the convection velocity is a function of the separation in time (or space). Convection 
velocities for the largest separations were estimated in this way at each r position 
and are summarized in Table 5.1. Shown are convection velocities normalized by the 
local mean velocity U{y ) and the average mean velocity U m = J^ 6 U(y ) dy. Note 

that at the locations nearer to the wall this convection velocity is generally closer to 
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TABLE 5.1 


Relative Convectiou Velocities at Large Separation 


y 

y + 

U c /U based 

on 

U c /U rn based on 

Reo 

Rrr 

R’ZZ 

Roo 

R rr 

Rzz 

0.813 

34 

0.78 

0.81 

0.89 

0.75 

0.78 

n iQ 

0.352 

117 

0.65 

0.54 

0.66 

0.74 

0.61 

C.'.. 

-0.352 

100 

0.79 

0.66 

0.79 

0.87 

0.73 

0.87 

-0.813 

29 

0.98 

0.97 

0.86 

0.85 

0.84 

0.74 


the local mean velocity than it is at locations farther from the wall. As explained 
by Sternberg (1967), large-scale structures extend over a significant portion of the 
flow; thus, their convection velocities are lower than the local mean Velocity in the 
regions away from the wall and higher near the wall. Consequently, there exists a 
point at which the convection velocity of large structures is equal to the local mean 
velocity. Apparently, the r locations closer to the wall are near this point because 
the ratio of the convection velocity to the local mean velocity is approaching 1. We 
neglected to calculate space-time correlations and convection velocities for y + < 30 
where we would expect the convection velocities to be greater than the local mean 
velocity. These arguments would also lead us to expect that the temporal spectra 
at the locations nearer the wall (Figure 5.40) would show better agreement with 
the streamwise spectra fer large wave numbers. As mentioned above, this appears 
to be the case. 


5.6 Turbulent Flow Structures 


In this section we will investigate the structure of an instantaneous velocity field 




from the present calculations. Presented arc contour plots of the velocity and vor- 
ticity fields, as well as plots of the velocity vectors. Note that in all the contour 
plots to follow, solid lines will depict a positive qu-.ntity and dashed lines a negative 
quantity. 

Figure 5.45 shows the streamwise velocity u' in (0, z)-planes near the concave 
wall at y + — 6.14 and near the convex wall at y + = 5.29. In these plots note the 
elongation of the positive and negative regions (streaks) in the streamwise direction 
in agreement with experimental observations (Kline et al. 1967). This is also 
consistent with the two-point correlation functions presented in §5.4. In contrast, 
the radial and spanwise velocities ( v ' and w') do not show this distinctive streakiness 
(Figures 5.46 and 5.47). The outstanding feature of the radial velocity is that 
strong fluctuations are concentrated in very small regions, with strong negative and 
positive fluctuations tending to lie next to each other in the spanwise direction. As 
will be seen below these regions contribute to the Reynolds shear stress. These 
characteristics are reflected in the large values of u-flatness F(v) at the walls, and 
the two-point correlation of v, R rr (6z) (see §5.1 and §5.4). The spanwise velocity 
also shows small regions of intense fluctuation, though not as small as the regions 
of intense radial fluctuations. These regions of intense w fluctuations tend to occur 
at the same positions as the regions of intense v fluctuations. Notice that near the 
convex wall, turbulence is less vigorous than on the concave side, as is evidenced by 
the concentration of contour lines. For each velocity component, the same contour 
levels are used for both sides of the channel. 

High and low speed streaks have been observed experimentally only near the 
walls. In agreement with this result, contour plots of u' at y = ±0.352 (Figure 
5.48) are not streaky. 

In order to determine the structures associated with the regions of intense v 
fluctuations mentioned above, we examine some enlargements of the area around 
one pair of intense regions near the concave wall (the framed area in Figures 5.45, 
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5.46 anil 5.47). Contours of u', v' and w' fluctuations in the ( 0 , 2 r)-planc at y f — 6.14 
are shown in Figure 5.49. The distance between the centers of the regions of positive 
and negative v' fluctuations (Figure 5.496) is 20 wall units, which is the approximate 
position of the negative minimum in R rr (5z) at this y location (see Figure 5.28). 
Comparison of the u' and v' contours shows that u' and v' in this region have the 
same sign and thus contribute positively to the Reynolds shear stress (recall that v 
is positive toward the wall on this side of the channel). Contours of u' and v' in an 
(r, $)-planc passing through this region (at the z location designated as A in Figure 
5.49) are shown in Figure 5.50. Here, a sharp gradient of u' is seen to occur along 
a front inclined upstream from the wall. A sharp gradient of this kind is t u e event 
detected by the VITA conditional sampling technique (Blackwclder & Kaplan 1979; 
Kim 1983). The inclined front observed here is in agreement with the observations 
of Praturi & Brodkey (1978) and the conditionally averaged velocity found by Kim 
(1984). In Figure 5.50c the velocity vectors projected into this plane are shown. 
In agreement with the observations of Praturi and Brodkey, a transverse vortex 
downstream of the front is barely discernible; this spanwise vortex is more apparent 
in other (r,0)-planes. 

In Figure 5.51, velocity vectors projected into several (r, z)-planes are shown; the 
0 locations of the planes are designated B through F in Figures 5.49 and 5.50. In 
plane B just downstream of the front, a single vortex can be seen very near the wall. 
The inflow and outflow sides of this vortex are seen in Figure 5.49 as the positive 
and negative regions cf v : fluctuation. This vortex is not discernible upstream oi 
this location. Note, however, that the region of negative v' in Figure 5.49 extends 
significantly upstream of the beginning of the vortex. Downstream, in plane ^ , the 
vortex is larger and farther from the wall. Farther downstream, in planes D through 
F, a second vortex of opposite sign appears and the pair lifts away from the wall. 
The vortex has moved beyond y + = 40 in a ."treamwise distance of approximately 
150 wall units. The fact that the vortex is much larger where it is distant from 
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Ihe wall is in agreement -with the behavior of the spnnwise two point correlation 
function R fr (bz) (see §5.4). Also, in §5.4 it was found that R r£ (6z) lid not indicate 
the presence of vortex pairs at y*" > 20. Note that the vortex pair in Figure 5.51 
would not contribute to the minimum in the spanwtac two-point correlation function, 
because the vortices are not at the same y location. Farther downstream, the vortex 
pair in plane F has lifted out of the domain, probably because of mutually induced 

velocity 

We emphasize that only a solitary vortex was observed near the wall in the (r, z)- 
p’ane passing through point B. Further downstream am- away from the wall, this 
vortex was joined by a vortex of opposite sign. In fact, in a random sampling cf 
(r, «)-planes, significantly more solitary vortices than vortex pairs were observed. 
In their dual-view, hydrogen-bubble flow-visualization, studies Smith & Schwartz 
(1983) observed twice as many patterns indicating the presence of streamwise rota- 
tion as those indicating streamwise counterrotating pairs. 

The streamwise and radial velocities u' and v' in the [9, r)-plane located at z = 
27r/3 (halfway across the domain in the z direction) are shown in Figure 5.52. In 
the v! plot, the most striking feature is that near the wall the structures (fronts) 
are inclined in the downstream direction, at a small angle from the wall. A similar 
behavior was observed by Moin & Kim (1982) in their calculations, by Rajagopalan 
& Antonia (1979) and Kreplin & Eckelmann (19796) in their experiments in a plane 
duct, and by Brown and Thomas (1977) in a boundary layer. The contours of v' 
do not exhibit this feature. Note that near both walls the regions of positive u' 
fluctuations correspond to fluid coming toward the wall ( v ' < 0 at the convex or 
lower wall, and v' > 0 at the concave wall). Also, the turbulence activity is much 
greater near the concave wall. 

In Figure 5.53 the streamwise velocity in the (r, z)-plane at 6 = 0.08 (halfway 
across the domain in the 6 direction) is plotted. Near both walls the :s.t®r mting 
regions of high and low speed fluid are apparent. Also, these regions a rs z afined 




to locations near the walls, where they can be counted to obtain the mean streak 
spacing. On the concave wall, the streak spacing is found in this way to be A ' — 107, 
and on the convex wall it is A f = 02 or 108, depending on how the streaks arc 
counted. These values are in very good agreement with those obtained from the 
two-point correlation function Poo(fiz), (see §5.4) and the generally accepted value 
of 100. Contours of the velocity components near the walls (y + < 50) in the (r, z)- 
plane have been examined and found to be in agreement with the results of Moin 
& Kim (1982). 

Streamwisc vorticity contours in the (r, ;z)-plane, both near the wall and across 
the entire channel, are shown in Figure 5.54. Note that the streamwisc vorticity 
is concentrated near the walls, as was indicated by the mis vorticity profiles (see 
§5.1), and that the concave side of the channel appears more turbulent than the 
convex side. Near the v we see that there are actually two regions in which 
the streamwise vorticity is concentrated. One is very near the wall (y 4 < 5), and 
the other a little away from the wall (c ntered about y £=• 20). These two regions 
correspond to the wall maximum and the local maximum which were observed in 
the rms streamwise vorticity profiles (see §5.1). Also notice that the regions of 
intense st’-eamwise vorticity at the wall tend to occur beneath intense regions away 
from the wall and are of opposite sign. This is the effect of the no-slip boundary 
condition on vortical structures slightly away from the wall. This change of sign is 
responsible for the minimum in the rms streamwise vorticity profile at y + — 5. 
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6. Summary and Conclusions 


There are specific conclusions to be drawn from this work in the areas of (*) the 
development of the numerical methods, (it) the structure of wall-bounded turbulent 
flows, and (tit) the effects of longitudinal curvature on turbulence. These will be 
covered in the following paragraphs. 

A spectral method for the solution of the Navier-Stokes equations for plane- 
channel flow and for flow between concentric cylinders has been developed. The 
method is based on divergence- free vector expansions and quasi- orthogonal func- 
tions, which offer advantages in memory requirements and computing speed. In 
addition, the method treats the boundary conditions and the continuity constraint 
exactly. This new method has been tested by computing the various states of Taylor- 
Couette flow. The results are in excellent agreement with available experimental, 
theoretical, and other numerical results. 

Mildly curved ( 6/R = 0.013), turbulent channel flow has been successfully sim- 
ulated, demonstrating the feasibility of computing fully developed wall-bounded 
turbulent flows without subgrid scale models. The computation does, however, ex- 
hibit marginal resolution, which is most apparent in the projection ~rror. This error 
is potentially significant near the walls, though its contribution to the Reynolds- 
stress balances is negligible and the details of the computed near-wall structure 
of the flow are in agreement with available experimental results. Further study is 
need. ' to determine what effect the projection error has on other properties of the 
flow ( e.g ., higher order statistical correlations), and what, if anything, car be done 
to alleviate it. 

A technique has been developed for calculating the temporal spectra and corre- 
lations with minimal cost in computation and computer memory. The technique 
relics on the smoothness of the energy spectrum in a turbulent flow to allow accurate 
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interpolatior of sparsely dist ributed computed spectral data. This is a promisiug 
method for obtaining temporal information from a turbulence simulation. 

In the numerical simulation of curved, turbulent channel flow, the overall agree- 
ment of the mean-velocity profile and higher order turbulence statistics with avail- 
able plane-channel data was good, except for the known effects of curvature. The 
present simulation ha? ; rovided a large amount of data concerning various features 
of wall-bounded turbulent Dows that are qualitatively unaffected by curvature. Near 
the walls, the flow was found to be dominated by alternating high and low speed 
streaks, elongated in the streaniwise direction with mean spanwise spacing of about 
100 wall units, in agreement with experimental observations. This value was ob- 
tained from both the instantaneous velocity fields and from two-point correlation 
functions. The velocity normal to the wall, however, is dominated by small spots of 
intense fluctuations, rather than streaks. Away from the walls, convection velocities 
of large-scale structures were found to be significantly lower than the local mean 
velocity, also in agreement with experiment.al observation. Althought streainwise 
vortices were observed in the vicinity of the wall, the two-point correlation functions 
showed little evidence of vortex pairs. Solitary vortices were observed to be more 
prominent near the wall than were vortex pairs. 

Many features of the present simulation are in excellent qualitative agreement 
with the large-eddy simulation of Moin & Kim (1982). For example, in agreement 
with their results, the streamwise rms vorticity profile near the wall exhibits two 
local maxima, one at the wall and the other slightly away from the wall. The max- 
imum away from the wall is attributed to the streamwise vortical structures which 
were observed in the vicinity of the wall. The wall maximum is due to the combined 
effects of the vorticity slightly away from the wall and the no-slip boundary condi- 
tions. Further, at + he wall the streamwise and spanwise vorticity are enhanced by 
the splatting effect through vortex stretching and the no-slip boundary conditions. 
Also in agreemer ; v ith Moin and Kim, but in disagreement with the experiments 
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of Kieplin & Eckelmann (1979a), the skewness of the velocity fluctuations normal 
to the wall, S(v), were found to become negative and to obtain a local ininimum 
slightly away from the wall. The terms in the balance equations for the elements 
of the Reynolds stress were in qualitative agreement with those of Moin & Kim 
(1982). In particular, it was found that viscous and turbulent d : ffusion acted to 
transport turbulence energy from the regions of maximum turbulence production, 
both toward and away from the wall. 

The effect of curvature on the mean velocity was found to be in agreement with 
experimental observations; namely, in local wall coordinates the mean- velocity pro- 
file sufficiently far from the wall ( yjr > 0.0015) on the concave side lies below the 
convex side profile and the plane-channel profile. By comparing the concave and 
convex sides of the channel, it was found that curvature had little effect on the 
turbulence statistics when they were normalized by local wall variables. The most 
notable exceptions to this were the Reynolds shear stress, the terms in the balance 
equation for the Reynolds shear stress, and the near-wall skewness and flatness fac- 
tors. On the concave side, the Reynolds shear stress is enhanced relative to the 
convex side shear stress. In the balance equation, the turbulent diffusion and veloc- 
ity pressure terms are significantly aitere"’. And, away from the wall, the so called 
“extra” term due to streamwisc turbulent diffusion is not negligible. 

A dominant feature of the computed flow was the Taylor- G dr ltler vortices. They 
made a significant contribution to the Reynolds shear stress on the concave side, 
and were found to affect the underlying turbulence by convecting it, and by inducing 
a strain field which locally affected turbulence production near the concave wall. 
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Appendix A 


Weak Formulation 


In this appendix wc consider the weak formulation of the forced Stokes equations, 

and its connection with the numerical method discussed in §2. We start with the 

classical form of the forced Stokes equations 

dv 
~dt 

V- v = °, ... ... (A-1) 


= -VP -VxVxv + f, 
v(« = 0) = v it 


v = 0 at the walls 
v periodic at periodic boundaries. 

The Reynolds number parameter is not necessary and is omitted for brevity. The 
viscous term, which is usually written V 2 , is written here as —V x V x , which is 
equivalent and more instructive in this case. Let u be any vector function such 
that V • u = 0 and u = 0 at the wall?; u must also be a function for which these 
conditions make sense, thus, 


ueV={u:u6 Hq, V • u = 0} . 


(A.2) 


When (A.l) is dot multiplied by u and integrated over the domain D, wc obtain 

/ d °-§f*=-/ 0 Vxu-Vxvdv + J u-fdv. (A. 3) 

where we have integrated the pressure term and the viscous term by parts. The 
pressvre term drops out. We define the inner product (v,u) to be J c v -u dv. The 
weak formulation is obtained by requiring v G *V and that (A. 3) be true for all 

u G V. 

v € V, v(t = 0) = v, , 

(u,^) = -(Vxu,Vxv) + (u,/) , (A.4) 

Vu e V • 
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If v is a solution to (A . 1) , then it must be a solution to the weak problem (A. 4). 
We will now show that it is the unique solution to (A. 4). Suppose there were a 
second solution to (A. 4), call it w. Then we would have 


( u ,^)=-(Vx u ,Vxw) + (u,f). 
Subtracting (A. 5) from (A. 4) and letting A = v - w we obtain 


(A.5) 


( u ,^) = -(Vx u ,Vxd). 


(A.6) 


At any time A(t) £ V, so 


Now 


and 


(A,~) = ~(VxA,VxA) . 

(A.7) 

(Vx^,Vx4)= [ ( V x A) 2 dv > 0 , 

(A.8) 

Jd 


a a , r a a i a r 9 , 


•m ) = J D A -m- dv= 29tJ D A dv - 

(A .9) 

§if D A * dv ~ 0 ’ 

(A.10) 


Thus, for all times, 


and since A = A 2 = 0 at t — 0 (v and w have the same initial conditions), 

f A 2 dv — 0 , (A. 11) 

Jd 

for all times. Therefore, since the square integral of A is zero, A is zero for all time 
and the solutions v and w are the same. 

We have shown that if there is a solution to (A.l) it is the unique solution to 
(A. 4), and is, therefore, also the unique solution to (A.l). However, it is possible 
that (A. 4) will have a solution when (A.l) has none. This can occur when the 
solution to (A. 4) does not have a second spatial derivative, so that the derivatives 
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in the viscous term of (A.l) are not defined. In this case, the solution of (A. 4) is 
often referred to as a weak solution to (A.l). Note also, that we need not have 
required that u G V. If we require (A. 2) to be true for all u in some space of which 
V is a subspace (as long as the integrals exist and the integrations by parts can be 
performed), the uniqueness result would still hold. 

We will assume that € ”V, which is an additional regularity condition. At any 
time we may consider (A l) as an equation lor given v and f; with the weak 
formulation, 


ov ,, 

*r eV ' 

dv 

(u, -gf) = ( u , — Vx V xv + /) , 

Vug!/. 


(A.12) 


This is a useful point of view for separating the spacial and temporal discretizations, 
by arguments similar to those used above, we find that (A.12) uniquely determines 
Note that the weak formulation requires that the pressure gradient term be 
orthogonal (».e., the inner product is zero) to all u € V; thus, ^ is the orthogonal 
projection onto ~\) of — VxVxv J-f. This is the projection mentioned in §2.1. There 
is the possibility that A.12 will not have a solution, which would occur if VxVxv+f 
were incompatible in the sense described in §2.2. 

Now consider the approximation to (A.12), 


^ i 

)=(u N .-VxVx y + /), 
Vu e V , 


(A.13) 


where V N is an iV-dimensional subspace of M . This is a special case of the weighted 
residual method presented in §2.1. Let e N be the error — §7^1 by letting u in 
(A.12) be and subtracting (A.13) from (A.12) we find that (u N ,e N ) = 0. 
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I 


Consider another approximate solution, %yf N G V N , with error e /N = e N +- 7 , 
where 7 G V N . We are interested in the Z «2 norm of 

(t'V") = («"«") + 2(^,7) + (7,7) • (A. 14) 

Since (c n , 7 ) = 0 and ( 7 , 7 ) > 0, we get 

(«",«")< (•'V") (A. 15) 

for all other approximations ■ Thus, the solution of A. 13 minimizes the Z >2 
norm of the error of the approximation to in "V . 


j L 
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Appendix B 


Saved Velocity Fields 


Velocity fields from the present calculation were saved every 200 time-steps (0.1 
time units) .This was accomplished by saving the restart files of each run. The 
first record of these files contains six integers that describe the data. They are 
T3TE P , NR = 96, NTHl = 128, NZl - 128, MR = 5, and MZ = 4; which are 
(in order), the time-step number, the number of planes in the r direction minus, 
one, the number of Fourier modes in the 9 direction, the number of modes in the 
z direction, the drawer size in the r direction, and the drawer size in the z direction 
(see §3 concerning drawers). Record two contains a single word used in computing 
the mean pressure gradient. The next NR/MR + 1 = 20 records contain the first 
column of drawers of the VDATA data base, which are followed by a single record 
containing a slab of the RDATA data base. This pattern of a column of VDATA 
drawers followed by an RDATA slab is repeated NZ1/MZ = 32 times. 

Each VDATA drawer is a complex array DATAV dimensioned as 
DATAV[NTKl/2,MZ,5,MR], which contains the Fourier coefficients of the veloc- 
ity field and its derivatives with respect to r. The first index indicates the 9 wave 
number of the coefficient (kg = 27r(Ii — 1)/0.16). The second index is for the z wave 
number, 


k t 


§ [h + (ZGROUP - l)MZ/2 - 1) h < MZ /2 

- §[I 2 - MZ/2 + (ZGROUP - l)MZ/2 - 1] I 2 > MZ/2 ’ 


where ZGROUP is the column number of the drawer. The third index determines 
the velocity component (1 — * vg, 2 — »v f) 3 — 4 — » 5 — ► )• And finally, 

the last index indicates the r plane (NRPLANE = I 4 + (RGROUP — 1)MR), where 
RGROUP is the row number of the drawer. The y location of the r plane is given 
by y = cos (x (NRPLANE - 1)/NR). Data in the VDATA drawers has been packed 
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to half precision and thus each record actually contains NTHl * MZ * 5 * NR/2 
physical words which must be unpacked before they are used. The slabs of the 
RDATA data base are useful only for restarting the computations, so these records 
should be skipped when using the data for another purpose. 

The restart files described above have been saved on XIOP tapes on the Cray- 
XMP at Ames Research Center using the STAGEX utility. The files are named 
‘CURVE002i’, where the “i” is either A or B. The individual restart files are dis- 
tinguished by their edition numbers. Table A.l lists the tape names and the files 
they contain. The owner of the tapes and the ID of the files is STTRDM. 
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Table A.l 


Restart Tapes 


Tapes 

Suffix x 

Edition numbers 

CURVE002ATAPE0M1 

A 

17, 18, 19, 20, 21 

CURVE002 ATAPEOO 1 

A 

22, 23, 24 

CURVE002ATAPE002 

A 

25, 26, 27 

CURVE002ATAPE003 

A 

28, 29, 30 

CURVE002ATAPE004 

A 

31, 32, 33 

CURVE002ATAPE005 

A 

34, 35, 36 

CURVE002ATAPE006 

A 

37 

EMERGENCY, ED=1 

A 

38, 39, 40, 41, 42, 43 

CURVE002ATAPE007 

A 

44, 45, 46 

CURVE002ATAPE008 

A 

47, 48 

CURVE002ATAPE009 

A 

49, 50, 51 

CURVE002ATAPE010 

A 

52, 53, 54 

CURVE002ATAPE011 

A 

55, 56, 57 

CURVE002ATAPE012 

A 

58, 59 

CURVE002ATAPE013 

B 

60, 61, 62, 63 

CURVE002ATAPE014 

B 

64, 65, 66, 67 

CURVE002ATAPE015 

B 

68, 69 

CURVE002Aj..iPE016 

B 

70, 71, 72, 73 

CURVE002ATAPE017 

B 

74, 75, 76 

CURVE002ATAPE018 

B 

77, 78, 79, 80, 81 
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Reynolds Stress Equations in Cylindrical Coordinates 


In this appendix the Reynolds-stress balance equations in cylindrical coordinates 
are presented. They arc derived from the balance equations in general tensor nota- 
tion, which may be obtained easily from the Navier-Stokes equations; they are 


du* u* 
dt 


= — U l j — (utnWj + u l u'U*^j — yu } u l u'^j ^ 

- (a ’ 1 (pu‘) [ + a ' 1 (pu 1 ) J + p + s’V,) 

+ (53) - («W + • 


(C.l) 


Here, U' is the mean velocity and it- 7 is the fluctuating velocity. The superscripts 
are contravariant indices, and the subscripts are covariant indices. Subscripts fol- 
lowing a comma indicate the covariant derivative. The summation con/ention is 
implied, and the overbars indicate ensemble averages; g t} is the contravariant met- 
ric tensor which in cylindrical coordinates is diagonal (since the coordinate system 
is orthogonal). For the cylindrical coordinate system with coordinates r, 0, and z , 
the contravariant elements of the metric tensor are given by 


9 = 


f; 


0 

l 

r* - 

0 



(C.2) 


Terms on the right-hand side of equation (C.l) are given the interpretations (in 
order of appearance) convection (referred to as Cii in the discussion below), pro- 
duction [Pij), turbulent diffusion ( TDij ), pressure diffusion ( ), pressure strain 
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correlation ( PSij ), viscous diffusion and viscous dissipation 

Thus, the transport equation for the Reynolds stresses ( R tJ ) can be written 


= C„ + K + TDa + PD„ + PS„ + ±-VD „ + ifl,, . 

(0.3) 

Each of these terms has been expanded in cylindrical coordinates and is presented 
in the equations that follow. In these equations, subscripts refer to the component, 
thus u r is the fluctuation part of the radial velocity and C r o is the convection term 
in the equation for u r uo (this is not to be confused with the general tensor notation 
used in (C.l)). A subscript following a comma indicates a derivative (not a covariant 
derivative); thus (u r u o) e ~ 9 • 

The convection terms C t] expand to 

Crr = -U r (u r U r ) f - — ((UrUr),* ~ 2u r U 0 ^ + Ux (u r U r ) f J 

(C.4a) 

Coe = -Ur {uou 0 ) <r - --- (jiueuo) g + 2u r tq>) - U x {u o u 0 ) t ; 

(C.4b) 

U* 

Cmx = —Ur {u z U,) (ti,ti,) 0 - Ux (llj Uj.) ; 

V 

(C.4c) 

CrO = -Ur {u r Uo) >r - — (^U r U 0 ) $ + U r U r ~ UflUs) ~ U M (u r Ue) 

, ; (C.4d) 

Cox = -Ur ( u 0 Ux) <r - — ({ueUx} <0 + U r u* ) - U, (u 9 u,) ( , ; 

(C.4e) 

C T x ~ ~U r {‘U r Ux)' T — —— (^{'dr u x) t o ~ u O u x^j ~ U z (lt r tij.) iX . 

(C.4f) 

The production terms expand to 


2 

P rr — —2Ur,rU r U r — ( Ur,0 ~ U$) V r Uo “ 2t/ r(> U r Ujr J 

(C.5a) 

2 

Poo = -2 Uo, r U r Uo (Uo ,0 + Ur) UoUQ - 2U 0i xUoUx ; 

(C.5b) 


H 

A 
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2 

P" ~ -2U 2 , r u~uZ — -Ug'OWu* - 2 U t ,*ujuz ; 

T 

PrO — ~ U T ,rU r Uo Ug >r U r U T — — {U r< $ — l Jo) UqUq — 

T 

- Ur, z uoul - U 0 , x ujujj ; 

Pot — Uo,r^r^z ~ — (Uo,0 ~F Ur) ti0U T — 

— Ue^u t u t U XiX t iou x ) 

Prz = U r . r U T U x — U XtT U r U r — — (U r ,6 — Uq) UoU z — 

r 

— i /r,iU,u t — U x<x u r u x . 


The turbulent diffusion terms TD,j expand to 

TD rr — (ru r u r ti r ) r - ^ ^(u r u,.ue) i() - 2u T u 0 Ug j 

TD 00 — - - (rU r UgUo) r — - ^(uoUpti'^) s - 1 - 2u r Uotig^ 

TD xx r= -i (ru r ii,u,) tr - £ (u7«7u7) >tf - (u^r,u7) j 

TD r0 ~ - (ru r u r u g ) r — ((ur 'o«„y * + n r u- u'g - 

- (u r u 0 u x ) z ; 

TD 0x — - ~ (rUrUgUi) r - - ^(u 0 U^u 7 ) g + UrlloU^j 

PD rz ' “ (fU r Urtt«) r — - l^(u r UgU x ) 6 • y lgUgU x ) 

The pressure diffusion terms PDij expand to 

rr = (fcp), r - w) ; 


(0.5c) 

; ( Uo,e + U r ) t bud 

(C.5d) 

Kr 

-UzoUgUo 

r 

(C.5e) 

;U x ,ou7uo 

(C.5f) 

I ~ ( u rU r u x ) ^ ; (C.6a) 

) - (u^u;) t ; (C.6b) 
r ; (C.6< ) 

UgUgUoj 

(CM) 

- (u 0 u x u x ) tM ; (C.6e) 

- (u t u 7u7) t . (C.6f) 


PD, 


(C.7a) 
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PDou ~ ~~ ((^«p),w + «rp) ; 

(C 7b) 


PD tz = -2{u~p) g ; 

(C.7c) 


PZ) r0 - -- (ru«p) r - - ((uTp) o - 2uiTp) ; 

(C.7d) 


PD 0 ; = - («op) iZ + £ («7p) )fl ; 

. (C.7e) 


PD ri = - - ((ru,p) r - u*p) - (u^) z . 

(C.7f) 

The pressure strain terms PS,j expand to 



PS rr = 2pu r>r ; 
2 

(C.8a) 


PS## = — p (ttfl.o + u r ) ; 

(C.8b) 


PS Z2 = 2 pu ZtZ ; 

(C.8c) 


PS r0 = P ^ * (llrf - Ug) + Ug^ ; 

(C.8d) 


P'S’ci = p > 

(C.8r) 


PStz — P ('Ur.i T . 

(C.8f) 

The viscous diffusion terms V Dij expand to 


VD rr = 

f ( r K“r),r) r + ^2 (i U »-«r) >e£ , ~ 4 (u r Ug) e + 2 (uyUj, - 

- U r ti r )j 


+ K«r) t „ ; 

(C.9a) 

V Dqo = 

- (r ( u o«fl), r ) r + ^2 ((«eWo),ee + 4 (“rM.o - 2 (u^u* - 

- u r u r )) 


+ ( u o«o) ; 

(C.9b) 


c-d- 
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GRlll ... 

OE FOGri i^^ r , 

VDzz = * (r(«r«:) <r ) r + (u7«7) >oe + (“!«;)„ ; 

VD r0 = * (r (t vuo) r ) ^ M 4 2 (u7u;) 9 - 2 («W) 9 

+ (uTSo),,- ; 

VD 0z = ^ (r (tt^tir) r ) 4- ^ ((tto«7), fia + 2 (tvul)* - SfltZl) 

+ ; 

VD rz = X - (r (u;tZ7) r ) ^ 4 — ((u^I 7) (M - 2 (u^),, - u^l) 

+ (tvu^),„ • 

The viscous dissipation terms Di } expand to 



(C.9c) 

4UrUo^ 

(C.9d) 

(C.9e) 

(C.9f) 

(C.lOa) 

(C.lOb) 

(C.lOc) 

(C.lOd) 

(C.lOe) 

(C.lOf) 
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Figure 1.1. Mean velocity profile for mildly curved channel How Re, = 1250. 

(Hunt and Joubert 1979). concave walll, convex wall, 

plane channel. 
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Figure 3.1. Flow of externally stored data in program CURVE. RDATA and 
VDATA are independent data bases, and PASSl and PASS2 are the 
two passes through the database each timestep (see §3.2). 
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Figure 4.1. Secondary flow stream-function contours for Taylor-Couette flow, 
(Re = 195, A = 2.004, r} = 0.95). Solid lines are positive contours; 
dashed lines are negative contours. 



Figure 4.2. Contours of axial velocity at r = 0.822r o for fully developed wavy 
vortices (Re = 458, A = 3.0, 77 — 0.868, m = 6). 
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Figure 4.3. Contours of axial velocity at r = 0.986r o for modulated wavy ^ortex 
flow (Re = 1300, A = 2.r ", m : = m 2 = 4, n = 0.877); (a) t = 12.59, 
(b) l = 11.58. 
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Figure 5.6. Rcyuolds stresses (a) in global coordinates (— uv), (— uo), 

viscous plus turbulent stress, equilibrium stress; (b) u'v' 

(c) utJ in local wall coordinates, concave side, convex side, 

o plane channel data from Eckelmann (1974). 
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Figure 5.7. Correlation coefficient with and without the contribution of the 
Taylor-Gortler vortices. ([wj 17 *)' 
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Figure 5.14. Secondary flow streamlines of the Taylor G or tier vortices, (a) not av- 
eraged over mirror image, (b) averaged over mirror image. 



Figure 5.15. Variation of the wall shear stress in the spanwise direction, con- 
cave wall. convex wall. 
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convex wall 


Figure 5.17. Contours of uv in the ( r,z ) j ie. Contour levels incremented by 0.1. 
Innermost positive contour is 1.55u^. 
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Figure 5.18. Contours of the diagonal elements of the Reynolds stress tensor due to 
the underlying turbulence in the (r, z)-plane. (a) ujj 2 — ttjp, (b) uj. 2 — u' 2 , 

(c) u'J 2 — ujp. Contour levels in (a) incremented by 0.2, in (b) and 
(c) by 0.05. Innermost negative contour in region “C”, (a) -1.9u?, 
(b) -,08u?, (c) — .22uJ. 



Figure 5.19. Contours of the Reynolds shear stress due to the underlying turbu- 
lence in the (r, z)-planc. Contour levels incremented by 0.05. Inner- 
most negative contour in Region “C” 0.18u^. 
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Figure 5.20. Terms in the balance of u 2 in local wall coordinates, (a) concave 
wall, (b) convex wall. Legend: production; radial tur- 
bulent diffusion; streamwise turbulent diffusion; velocity- 

pressure gradient; dissipation; radial viscous diffusion. 
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Figure 5.26. Terms in the balance of 2 q 2 = u 2 -f- u 2 + in 2 in local wall coordinates, 
(a) concave wall, (b) convex wall. See caption of Figure 5.19 for 
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Figure 5.31. Isocorrelation contours in the (#,;:)■ plane at y = .966, y h — 6.13, 
(a) Root (b) R T t, (c) R Z z • Contour levels are incremented by 0.1. 
Domain is 720 local wall units in 0 and 170 units in z. 
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Figure 5.32. Isocorrclation contours in the (0, z)-piane at y = .352, y + = 117, 
(a) Roo, (b) Rrrt ( c ) Rzz • Contour levels incremented by 0.1. Domain 
is 670 global wall units in 0 and 440 units in z. 










Figure 5.33. Isocorrelation contours in the (0,z)-plane at y = -.966, y + = 5.28, 
(a) Rflo, (b) R rr , (c) R zz . Contour levels incremented by 0.1. Domain 
is 605 local wall units in 0 and 405 units in z. 










Figure 5.34. Streamwise one dimensional energy spectra. 
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Figure 5.36. See next page for caption 



















Figure 5.41. Iso correlation contours in the (0,t)-plane at y = .813, y + = 34 
(a) Roo, (b) R rr , (c) R ig . Contour levels incremented by 0.1. Do 
main is 0.5 6/u T in time and 12.85 in 6. 








Figure 5.43. Isocorrelation contours in the (0, £)-plane at y = -.352, y + = 100, 
(a) Boo , (b) Rrr , (c) R z: . Contour levels incremented by 0.1. Domain 
is 0.56/u r in time and 12.68 in 6. 
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Figure 5.45. Contours of streamwise velocity tt' in the (0, z)-plane, (a) near the 
concave wall, y+ = 6.14, (b) near the convex wall, y + = 5.29. 
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Figure 5.47. Contours of spanwise velocity w' in the (0,z)-plane, (a) near the con 
cave wall, y f — 6.14, (b) near the convex wall, j/ H — 5.29. 







(a) 


ORIGINAL FAOiE tSI 
OF POOR QUALITY 



(b) 




Figure 5.48. Contours of streamwise velocity u' in the (0, z)-plane, (a) near the 
concave wall, y = .352, j/ + = 117, (b) near the convex wall, y = —.352, 
y + = 100. 
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Figure 5.49. Contours of (a) u', (b) v' and (c) w' in the (0,z)-plane near the concave 
wall, y + = 6.13. Enlargement of the framed region in Figures 5.45a, 
5.46a and 5.47a. The domain is 135 wall units in the z direction and 
404 wall units in the 9 direction. 
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(b) concave wall 



(c) concave wall 



Figure 5.50. Contours of (a) u' , (b) v' , in the (r,0)-plane and (c) velocity vectors 
projected into the (r, 0)-plane at the z location marked “A” in Figure 
5.50. The domain is 404 wall units in the 0 direction, and extends to 
— 52.7. 
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Figure 5.51. Velocity vectors projected into (r, z)-planes at the 6 locations marked 
“B” through “F” in Figures 5.49 and 5.50. The domain is 135 wall 
units in the z direction. 
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Figure 5.52. Contours of (a) u' and (b) v' in the (0,r)-plane at z = 2n/3. The 
radial direction has been magnified by a factor of 2.5. 
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Figure 5.54. Contours of streamwise vorticity u>' 0 in the (r, e)-plane at Q — 0.08; 

(a) across the entire channel — 1 < y < 1, (b) near the concave wall 
0 < y + < 49, (c) near the convex wall 0 < y J < 42. In (b) and (c) 
the radial direction has been significantly magnified. 









